•flrersa.  ^1(2 -  q  ~7  -cno? 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704  0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing 
the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  forTsducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information 
Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1 204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank) 


2.  REPORT  DATE 

August  1996 


4.  TITLE  AND  SUBTITLE 

Active  Control  of  Aerolasticity  and  Internal  Flows  in  Turbomachinery 


3.  REPORT  TYPE  AND  DATES  COVERED 

Final  Technical  Report  1  Nov  92  to  30  Apr  96 


5.  FUNDING  NUMBERS 

F49620-93 - 1 -0032 


6.  AUTHOR(S) 

Alan  H.  Epstein 
James  D.  Paduano 
Edward  M.  Greitzer 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Gas  Turbine  Laboratory 
Department  of  Aeronautics  and  Astronautics 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139-4307 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
AFOSR/NA 

110  Duncan  Avenue,  Rm  B115 
Bolling  AFB,  DC  20332-8050 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

F49620-93- 1-0032 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


13.  ABSTRACT  (Maximum  200  words) 

The  research  conducted  was  focused  in  two  areas.  The  first  is  active  control  of  rotating  stall  when  inlet  distortion  is  present 
This  research  was  carried  out  on  the  low-speed  3-stage  active  control  research  compressor  at  MIT=GTL.  The  second  area 
of  research  was  active  control  of  surge  in  centrifugal  and  axi-centrifugal  engines.  Stabilization  of  such  engines  beyond  their 
normal  stability  boundary  requires  control  of  ID  oscillations  -  which  lead  to  surge  -  rather  than  the  higher-dimensional 
oscillations  which  lead  to  rotating  stall.  Typically  surge  control  alone  allows  uncontrolled  rotating  stall  modes  to  go 
unstable,  which  would  be  debilitating  (and  probably  eventually  lead  to  surge)  in  an  axial  compression  system.  Rotating  stall 
in  centrifugal  compressors,  however,  is  quite  different.  If  one  can  prevent  surge  in  a  centrifugal  machine,  rotating  stall 
causing  gradual,  recoverable  performance  degradation  of  the  compressor.  Thus  active  control  of  surge  alone  can  increase 
engine  operating  range  in  engines  containing  centrifugal  compressors. 


19971217  045 


14.  SUBJECT  TERMS 

IS.  NUMBER  OF  PAGES 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

18.  SECURITY  CLASSIFICATION 

19.  SECURITY  CLASSIFICATION 

20.  LIMITATION  OF 

OF  REPORT 

OF  THIS  PAGE 

OF  ABSTRACT 

ABSTRACT 

UNCLASSIFIED 

UNCLASSIFIED 

UNCLASSIFIED 

UL 

rmm  quale??  msmurm  & 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/D10R,  Oct  94 


GAS  TURBINE  LABORATORY. 
DEPARTMENT  OF  AERONAUTICS  AND  ASTRONAUTICS 
MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
CAMBRIDGE ,  MA  02139-4307 


Abstracts  for  Grant  #F49620-93-l-0032 
entitled 


ACTIVE  CONTROL  OF  AEROL ASTICITY  AND  INTERNAL 
FLOWS  IN  TURBOMACHINERY 


prepared  for 
AFOSR/NA 

110  Duncan  Avenue,  Suite  B115 
Bolling  AFB,  DC  20332-0001 

Dr.  James  M.  McMichael,  Program  Manager 


PRINCIPAL 

INVESTIGATOR:  Alan  H.  Epstein 

Professor  and  Director  of  Gas  Turbine  Laboratory 
Department  of  Aeronautics  &  Astronautics 


CO-INVESTIGATORS:  James  D.  Paduano 

Associate  Professor 

Department  of  Aeronautics  &  Astronautics 
Edward  M.  Greitzer 

H.  N.  Slater  Professor  and  Associate  Head 
Department  of  Aeronautics  &  Astronautics 

August  1996 


aCTTVE  CONTROL  OF  AEROIASTICITY  AND  INTERNAL  FLOWS  IN  TURBOMACHINERY 
ACnVE  CC^dn  RNL  Gr^J.  D. Padu.no,  C M.  VanSdaaUcwyk, F.'Alessa,  B.  Com 

MIT  Gas  Turbine  Laboratory,  Cambridge,  MA 

_  ,  ,  j  <  •„„  1  qq<  i  qt)6  was  focused  in  two  areas.  The  first  is  active  control  of  rotating  stall 
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Table  I:  Stable  Operating  Range  Increase  of  Various  Control  Laws,  0.8  Dynamic  Head  Distortion 


Type  of  Control 

Range  Increase 

Linear  Quadratic  Gaussian 

1.5% 

Cross-Coupled  Harmonic  Feedback 

3.0% 

SISO  Harmonic  Feedback 

2.2% 

Distributed  Feedback 

3.7% 

Two  other  control  laws  were  considered,  which  are  independent  of  distortion  location  and  magnitude.  These  are 
described  here,  and  compared  in  Table  I. 

The  first  ‘insensitive’  control  law  is  simply  a  SISO  controller  -  each  measured  SFC  is  fed  back  to  a 
corresponding  SFC  of  actuator  deflection.  This  allows  one  to  separately  stabilize  each  of  the  unstable  modes  during 
clean  inlet  flow.  During  distorted  flow,  however,  SFC  coupling  degrades  the  performance  of  such  a  control  law.  This 
controller  is  at  one  end  of  the  spectrum  of  trading  off  performance  (range  extension)  for  robustness  (insensitivity  to 
distortion  location).  A  simple  method  to  improve  performance  is  to  experimentally  introduce  and  tune  cross-coupling 

terms;  this  however  makes  the  control  law  sensitive  to  the  location  of  the  distortion  pattern.  ° 

The  second  ‘insensitive’  control  law  developed  is  termed  distributed feedback.  This  control  law  measures  the 
entire  measured  shape  of  the  disturbance,  and  feeds  back  a  shifted  and  amplified  version  of  it  to  the  actuators.  This 
control  law  has  several  advantages:  It  is  a  fixed  rather  than  dynamic  controller,  thus  very  simply  implemented.  It  is 
insensitive  to  the  location  of  the  distortion.  It  simultaneously  considers  all  of  the  Fourier  harmonics,  and  utilizes  the 
measured  relative  phases  of  the  SFCs  to  direedy  determine  the  phase  of  the  actuation  -  thus,  according  to  simulation 
studies,  it  is  also  insensitive  to  the  type  of  distortion  introduced.  Finally,  this  control  law  performed  as  well  or  better 
than  any  other  control  law  tested.  Interestingly,  the  trade-off  between  range  extension  and  insensitivity  is  broken  by  the 
distributed  feedback  controller  -  it  is  both  the  best  performer  and  the  most  insensitive  to  distortion  type  and  location. 


2.0  Surge  Control  in  Centrifugal  and  Axi -centrifugal  Engines 

The  second  area  of  research  we  will  report  is  active  control  of  surge  in  engines  with  centrifugal  or  axi- 
centrifugal  compressors.  Stabilization  of  such  engines  beyond  their  normal  stability  boundary  requires  control  of  ID 
oscillations  -  which  lead  to  surge  -  rather  than  the  higher-dimensional  oscillations  which  lead  to  rotating  stall. 
Typically  surge  control  alone  allows  uncontrolled  rotating  stall  modes  to  go  unstable,  which  would  be  debilitating 
(and  probably  eventually  lead  to  surge)  in  an  axial  compression  system.  Rotating  stall  in  centrifugal  compressors, 
however,  is  quite  different.  If  one  can  prevent  surge  in  a  centrifugal  machine,  rotating  stall  causes  gradual,  recoverable 
performance  degradation  of  the  compressor.  Thus  active  control  of  surge  alone  can  increase  engine  operating  range  in 
engines  containing  centrifugal  compressors. 

2.1  Experimental  Setup  and  Model  Development 

The  Lycoming  LTS-101  engine  was  outfitted  with  an  improved  valve  which  injects  flow  near  the  leading  edge 
of  the  vaned  diffusers.  This  was  the  mechanism  used  previously,  but  the  current  setup  has  been  modified  with  a  larger 
area  valve  and  enlarged  holes  leading  from  the  annular  plenum  which  feeds  into  the  hollow  diffuser  vanes  (these  in  turn 
have  slots  which  lead  into  the  diffuser  passages).  These  modifications  allow  up  to  4%  of  the  compressor  mass  flow  to  be 
modulated  by  the  actuator,  which  has  a  bandwidth  of  approximately  350  Hz. 

The  injectors  modify  the  steady-state  operation  of  the  engine  because  there  is  a  desire  to  modulate  actuation  both 
positively  and  negatively.  Thus  a  mean  injection  rate  of  2%  is  used,  and  the  modified  engine  operation  is  studied. 
Specifically,  system  identification  of  the  surge  dynamics  is  conducted  at  various  mass  flows,  to  determine  whether  the 
expected  linearized  surge  dynamics  exist  and  lead  to  instability  of  the  system. 

Figure  5  shows  one  such  identification  result.  Also  shown  is  a  transfer  function  fit  to  the  data,  which  required  3 
poles  (eigenvalues)  and  3  zeros.  The  eigenvalue  at  30  Hz  is  predicted  by  our  surge  models  with  appropriate  model 
parameters.  The  other  dynamics  are  interpreted  as  longitudinal  wave  dynamics  (acoustic  modes)  which  couple  with  the 
compressor  djmamics.  Dynamics  such  as  these  are  also  predicted  to  occur,  by  more  sophisticated  modeling  methods. 

To  verify  that  surge  is  in  fra  driven  by  instability  of  the  eigenvalue  at  30  Hz  as  prediaed  by  the  model,  the 
engine  operating  condition  was  driven  down  to  the  surge  point  using  its  movable  nozzle.  Figure  6  shows  very  clearly 
the  three  resonant  modes  exhibited  by  the  system.  This  data  was  taken  at  a  corrected  mass  flow  very  near  the  surge 
point.  Surge  inception  itself  appears  to  involve  energy  being  fed  from  the  30  Hz  mode  to  a  23  Hz  limit  cycle  which  is 
very  small  amplitude  compared  to  surge.  This  frequency  eventually  develops  into  a  deep  surge. 

2.2  Control 

Control  law  development  is  currently  underway  based  on  the  models  described  above.  Two  types  of  control  law 
have  been  designed.  Linear  controllers  which  rely  on  the  identified  dynamics  are  currently  being  tested,  and  closed 
loop  identification  results  will  be  used  to  improve  these  designs.  Nonlinear  sliding  mode  controllers  have  also  been 
designed,  and  will  be  implemented  to  assess  the  additional  effeaiveness  available  from  using  larger  amplitude 
nonlinear  actuator  commands  to  stabilize  the  system  in  the  free  of  modeling  errors  and  disturbances. 


Figure  1  -  Experimental  Setup  for  Active  Control  of  Rotating  Stall  with  Distortion. 
Distortion  ring  is  on  the  left,  shaded  parts  can  be  rotated. 


Figure  2  -  Steady-state  compressor  performance  with  and  without  1.9  dynamic  head  distortion. 
Extended  Hynes-Greitzer  prediction  based  on  uniform  flow  curve  fit  and  geometry  only. 
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Figure  3  -  static  pressure  upstream  of  the  compressor  due  to  a  1.9  dynamic  head  distortion. 
Bias  of  ps=.03  between  predictions  and  measurements  has  been  removed.  The  distortion 
screen  blocked  the  flow  between  120  deg  and  240  deg,  shown  by  dotted  lines. 
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Figure  4  -  Sample  transfer  function  identification  and  prediction  comparisons,  between  IGV 
deflection  and  axial  flow  perturbation.  Phase  and  magnitude  were  compared  for  all  transfer 
functions  in  the  3-by-3  matrix  between  first,  second,  and  third  harmonic  inputs  and  outputs. 


Figure  5  -  Measured  engine  transfer  function  with  new  injection  actuation,  between  mass 
flow  injection  and  upstream  static  pressure.  Solid  line  is  a  fit  based  on  a  3  pole  3  zero 
transfer  function.  Pole  locations  were  verified  over  several  operation  points. 


Figure  6  *  Spectrum  of  upstream  static  pressure  perturbations  immediately  prior  to  surge. 
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Active  Control  of  Rotating  Stall  with  Inlet  Distortion 

by 


Christiaan  Mauritz  Van  Schalkwyk 


Abstract 

This  thesis  describes  the  first  use  of  active  control  to  extend  the  stable  operating  range 
of  a  multi-stage  axial  compressor  in  the  presence  of  circumferential  inlet  total  pressure 
distortion.  Three  control  strategies  at  different  levels  of  sophistication  were  examined: 

(i)  Constant  gain  control  laws,  experimentally  optimized  to  obtain  the  lowest  stalling 
flow  coefficients  for  circumferentially  uniform  flow,  were  found  to  be  effective  in  the  presence 
of  distortion.  Taking  into  account  the  coupling  between  harmonics  (which  occurs  as  a  result 
of  distortion)  further  increased  the  stable  operating  range.  The  range  extension  was  3%  for  a 
0.8  dynamic  head  distortion,  (ii)  A  new  constant  gain  control  law  with  a  single  spatial  phase 
shift  was  also  examined.  This  stabilized  the  compression  system  at  lower  flow  coefficients 
for  distorted  flow,  increasing  the  operating  range  by  3.7%.  (iii)  A  theoretical  model  of 
the  unsteady  flow  response  of  an  axial  compressor  to  inlet  distortion,  developed  by  Hynes 
and  Greitzer,  was  used  to  design  linear  quadratic  Gaussian  controllers.  These  model  based 
controllers  achieved  range  extensions  of  1.5%  and  1.1%  for  inlet  circumferential  distortions 
of  0.8  and  1.9  dynamic  head  respectively. 

System  identification  procedures  were  used  to  examine  the  dynamic  behavior  of  the 
disturbance  modes  in  the  compressor.  The  dynamics  of  small  velocity  perturbations,  mea¬ 
sured  here  for  the  first  time  in  nonuniform  flow,  were  accurately  predicted  by  the  model, 
extended  to  include  unsteady  viscous  effects.  Experimentally  measured  multi- input  multi¬ 
output  transfer  functions  confirmed  the  strong  coupling  between  harmonics  of  small  velocity 
perturbations,  also  as  predicted  by  the  model. 
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Notation 


Overbars  and  Tildes 


Overbars  denote  averaged  quantities,  spatial  or  temporal.  Tildes  are  used  to  indicate 
the  (real  or  complex)  spatial  Fourier  coefficients  of  the  corresponding  variable;  see  Equa¬ 
tion  (3.1)  p.48. 


Subscripts 


1  upstream 
4  downstream 
a  AGV 
c  compressor 

cn  nth  cosine  Fourier  coefficicient,  n  =  1, 2, . . . 
i  inlet,  isentropic,  imaginary, 
m  axial  measurement  location 
n  harmonic  number 
r  rotor,  real 
s  stator,  sound 
s  steady 

sn  nth  sine  Fourier  coefficicient,  n  =  1, 2, . . . 
t  total  or  throttle 
w  wheel 
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Superscripts 

(•)*  complex  conjugate 
(•)T  transpose 

(■)ff  complex  conjugate  transpose 

1  far  upstream  (see  Figure  3.1  p.51  for  station  numbers) 

2  AGV  inlet 

3  AGV  exit  =  compressor  inlet 

4  compressor  exit 

5  plenum 

6  throttle  exit 

Vector  notation 

Vectors  will  be  written  as  a  column  and  will  be  denoted  by  bold  lowercase  characters.  The 
n  x  1  vector  is  thus  written  as 

Xl 

xn  _ 

=  [xi ,  .  •  ■  ,  Xn]T 

Matrix  notation 

Uppercase  letters  will  be  used  to  denote  matrices,  the  corresponding  lowercase  letters  with 
subscripts  ij  will  be  used  to  denote  the  (i,j)  entry,  so  the  m  x  n  matrix  is  written  as 

on  aln 

A  =  •  •  (0.3) 

flml  ’  * 4  flmn 

Columns  of  the  matrix  will  be  denoted  by  the  vectors  oi,  ...  ,  an.  The  transpose  of.  the 
matrix  will  be  written  as  AT  and  the  complex  conjugate  (or  Hermitian)  transpose  as  AH . 


(0.1) 

(0.2) 
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Diagonal  and  block-diagonal  matrices  will  be  indicated  by 

A  =  dia.g[Ai,A2,...,An] 

where  A4  is  the  ith  element  on  the  diagonal. 
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ba  Equation  (3.37)  p.57 

Dg  Fourier  spatial  derivative  matrix,  Equation  (3.6)  p.50 
EXm  Equation  (3.45)  p.59 

AGV  flow  coefficient  relations,  Equation  (3.10)  p.53 
Equations  (3.19)  and  (3.20)  p.54 

F(-)  Fourier  convolution  matrix,  Equation  (3.3)  p.49  and  Appendix  A 
G{ioj)  compression  system  transfer  function 
i  y/—l 

k  controller  gain  for  distributed  feedback  controller 
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p  pressure 


Table  continues  on  next  page. 
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p^Pj  AGV  pressure  relations,  Equation  (3.12)  p.53 
Equations  (3.21)  and  (3.22)  p.54 
r  reaction 
f  mean  rotor  radius 
t  time 


it  velocity 

vi,V2,--.  eigenvectors  of  A  corresponding  to  eigenvalues  A*,  A2,  — 

V  matrix  of  eigenvectors,  Equation  (4.4) 
x  nondimensional  axial  location 
xm  nondimensional  upstream  sensor  location  (m=measure) 

6x  state  vector,  Equation  (3.42)  p.58  and  Equation  (3.68)  p.65 


P 

Pn  >  0nn 
6 

7 

Ai,  A2, . . . 
A 
P 

4> 

ip 

p 

e 

T 

Tf 

Tr,T3 

OJ 

Uc 

U)T 


spatial  phase  shift  of  distributed  controller 
spatial  phase  shift  for  nth  harmonic 
perturbation 
AGV  deflection 

eigenvalues  of  A  corresponding  to  eigenvectors  i>i,  i>2>  •  •  •  • 
matrix  of  eigenvalues,  Equation  (4.5) 
fluid  inertia,  see  p.  17  for  meaning  of  subscripts 
flow  coefficient 

compressor  inlet-total  to  exit-static  pressure  rise 

density  of  air 

angle  around  annulus 

nondimensionalized  time 

constant,  Equation  (3.59)  p.63 

time  constants  for  loss  dynamics,  Equations  (3.61)  and  (3.62)  p.63 
nondimensionalized  frequency 
cross-over  frequency 
rotor  frequency 
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Abbreviations 


AGV  actuator  guide  vanes 
DF  distributed  feedback 
DFT  discrete  Fourier  transform 
FIR  finite  impulse  response 
HF  harmonic  feedback 
HFC  harmonic  feedback  with  cross  coupling 
IMS  integrated  mean  slope 
CIMS  corrected  integrated  mean  slope 
LQG  linear  quadratic  Gaussian 
MIMO  multi-input  multi-output 

MIT  Massachusetts  Institute  of  Technology 
PSD  power  density  spectrum 
SISO  single-input  single-output 
SM  surge  margin 
SNR  signal  to  noise  ratio 
STD  standard  deviation 


XHG  extended  Hynes-Greitzer  model 
XHG-Euler  XHG  with  nonlinear  description  for  upstream  flow  field 


23 


24 


1  Introduction 


Compression  systems  axe  subject  to  instabilities  that  limit  the  range  over  which  safe  oper¬ 
ation  is  possible.  As  the  mass  flow  through  the  compressor  is  decreased,  the  peak  pressure 
rise  across  the  compressor  increases  until  a  point  is  reached  where  the  flow  through  the 
compressor  becomes  unstable.  Loss  of  stability  is  undesirable,  as  the  amplitudes  of  the 
unstable  oscillations  are  often  very  large  and  can  cause  damage  to  the  engine.  In  addition, 
the  loss  of  stability  is  accompanied  by  a  significant  loss  in  pressure  rise. 

Two  types  of  instabilities  are  commonly  observed  in  compression  systems.  The  first, 
called  surge,  is  a  system-type  instability  that  is  mainly  axisymmetric.  Sometimes  the  surge 
oscillations  axe  so  severe  that  flow  through  the  compressor  reverses.  The  second,  called 
rotating  stall,  has  the  form  of  a  stall  cell  that  rotates  around  the  annulus  at  a  fraction 
of  the  rotor  speed.  Often  rotating  stall  starts  out  as  a  small  perturbation  wave  travelling 
around  the  annulus,  which  grows  in  amplitude,  and  finally  develops  into  a  large  amplitude 
limit  cycle,  called  rotating  stall.  Due  to  hysteresis  the  only  way  to  stop  this  limit  cycle 
oscillation  is  to  increase  the  mass  flow  significantly  beyond  the  point  where  stall  occurred. 
Whether  the  compressor  surges  or  stalls  depends  on  the  specific  system  configuration  and  is 
determined  by  the  ratio  of  the  plenum  compliance  to  duct  inertia  as  shown  by  Greitzer  [17]. 
A  review  of  axial  compressor  stall  phenomena  is  given  by  Greitzer  [16]. 

Because  of  these  aerodynamic  instabilities  operation  near  the  unstable  region  must  be 
avoided.  However,  this  margin  of  safety  forces  operation  at  mass  flows  that  deliver  lower 
than  the  maximum  peak  pressure  rise.  The  surge  margin  (SM)  adopted  by  Aerospace 
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Recommended  Practice  committee  [42]  is  defined  by 

SM  =  {Rs  -  Rq)/Ro  (1-1) 

where  Rs  is  the  surge  pressure  ratio  and  Ro  is  the  operating  pressure  ratio.  A  typical  value 
of  SM  is  in  the  range  10  -  20%,  and  30  -  50%  of  this  number  is  accounted  for  by  nonuniform 
inlet  conditions  [50]  which  we  discuss  next. 

If  the  flow  through  the  compressor  is  nonuniform  the  pressure  rise  across  the  compressor 
decreases.  In  addition,  the  compression  system  goes  unstable  at  higher  mass  flows.  This 
nonuniform  flow  is  also  called  distortion. 

There  axe  several  causes  of  distortion.  During  takeoff  the  flow  enters  the  inlet  at  an 
angle,  resulting  in  nonuniform  total  pressure  around  the  annulus.  If  an  aircraft  is  flying 
at  high  angles-of-attack  the  flow  can  separate  from  the  inlet  lip  and  this,  too,  will  create 
nonuniform  flow  through  the  compressor.  Nonuniform  temperature  profiles  can  exist  in  the 
inlet  ducts  in  very  short  takeoff  and  landing  aircraft  when  hot  exhaust  gases  are  ingested. 
Whatever  the  cause,  distortion  adversely  affects  compressor  operation. 

The  distorted  flow  can  be  circumferentially  or  radially  nonuniform.  Reid  [41]  found  that 
circumferentially  nonuniform  flow  is  more  severe  so  we  will  concentrate  on  it.  Reid  found 
that  if  the  angular  width  of  low  inlet  total  pressure  is  increased  from  zero  there  is  an  initial 
steep  drop  in  the  exit  static  pressure  of  the  compressor  until  an  angle  is  reached  beyond 
which  there  is  little  change  in  the  exit  static  pressure.  This  angle  is  called  the  critical  sector 
angle  and  is  shown  in  the  graph  at  the  top  of  Figure  1.1,  taken  from  Reid  [41].  Reid  also 
found  that  if  the  distorted  region  is  split  into  several  smaller  regions  distributed  around 
the  annulus  but  the  cumulative  width  is  kept  constant,  the  biggest  loss  in  pressure  rise 
occurs  when  there  is  only  one  continuous  distorted  region.  This  is  shown  in  the  graph  at 
the  bottom  of  Figure  1.1. 


26 


(a)  Effect  of  Spoiled  Sector  Width 


®(2>  ® 


(b)  Contiguous  Spoiled  Sector  Width  is  Important 

Figure  1.1:  Effect  of  distortion  on  exit  static  pressure.  Taken  from  Reid  [41]. 


The  loss  in  pressure  rise  can  be  explained  by  a  very  simple  model.  Assume  that 
two  halves  of  the  compressor  annulus  are  operating  at  two  different  pressures  as  shown  in 
Figure  1.2.  In  this  figure  we  see  that  the  mean  peak  total  to  static  pressure  rise  Vts  is 
lower  than  that  of  the  compressor  operating  at  the  same  mean  mass  flow.  This  model, 
proposed  by  Mazzawy  [32],  is  often  called  the  parallel  compressor  model.  Although  the 
parallel  compressor  model  explains  the  loss  in  peak  pressure  rise  it  does  not  capture  all  the 
dynamics  in  a  compression  system.  Hynes  and  Greitzer  [22]  realized  that  to  take  distortion 
into  account  the  nonuniformity  must  enter  the  dynamics  in  a  nonlinear  way  and  developed  a 
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Low  Total 
Pressure 


High  Total 
Pressure 

Annular  Cross-Section 
Showing  180°  Distortion 


Figure  1.2:  Parallel  Compressor  Model.  Taken  from  [29]. 


model  to  account  for  distorted  flow.  Studies  by  Hynes  and  Greitzer  [22]  found  that  the  model 
correctly  predicts  the  trends  that  were  observed  by  Reid  [41].  Hynes  and  Greitzer  [22]  and 
Chue  et  al.  [5]  also  found  that  the  mass  flow  at  which  the  circumferentially  integrated  mean 
slope  (IMS)  of  the  compressor  characteristic  is  equal  to  zero  approximately  corresponds  to 
the  point  where  the  compression  system  looses  stability.  The  model  by  Hynes  and  Greitzer 
is  an  extension  to  the  distorted  flow  case  of  an  earlier  model  by  Moore  and  Greitzer  [34]. 
A  review  of  different  models  for  the  analysis  of  compression  system  instabilities  is  given  by 
Longley  [28]. 

The  model  by  Moore  and  Greitzer  [34]  predicts  that  a  small  amplitude  perturbation 
should  propagate  around  the  annulus  prior  to  stall.  McDougall  [33]  observed  these  waves  in 
a  low  speed  single-stage  axial  compressor,  and  Garnier  [12]  identified  them  in  both  a  single- 
and  a  three-stage  compressor.  Experiments  by  Day  [7]  showed  that  some  compressors  do 
not  exhibit  travelling  waves  prior  to  stall;  instead  a  short  length  scale  perturbation  rotated 
around  the  annulus  at  approximately  70%  of  the  rotor  speed.  The  reason  for  this  behavior 
is  not  known  and  is  currently  being  researched.  Day  and  Freeman  [8]  and  Tryfonidis  et 
al.  [46]  observed  that  travelling  waves  are  also  present  in  high  speed  compressor  prior  to 
stall.  These  are  encouraging  results  as  it  allows  experiments  to  be  conducted  on  low  speed 
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compressors. 


The  Hynes-Greitzer  model  predicts  that,  in  the  presence  of  distortion,  small  amplitude 
waves  will  travel  around  the  annulus.  But,  unlike  uniform  flow,  the  waves  are  not  purely 
sinusoidal  but  change  shape  as  they  propagate  around  the  annulus.  Longley  [25]  performed 
distorted  flow  experiments  and  observed  the  predicted  circumferential  variation  in  mag¬ 
nitude  of  the  propagating  disturbances.  He  also  derived  an  improved  stability  criterion, 
called  the  corrected  integrated  mean  slope  (CIMS)  criterion,  and  showed  that  the  ratio  of 
the  plenum  compliance  to  duct  inertia  derived  by  Greitzer  [17]  determines  the  applicability 
of  the  CIMS.  Longley  also  found  that  the  pressure  rise  predicted  by  the  model  is  lower  than 
the  experimentally  observed  pressure  rise,  and  attributed  the  discrepancy  to  the  large  inter¬ 
blade  row  gaps  and  sensitivity  to  inlet  swirl  of  the  experimental  compressor.  A  method  to 
determine  the  rotor  fluid  inertia,  an  important  parameter  in  the  model,  from  steady  state 
velocity  profiles  was  also  given  by  Longley  [25]. 

To  reduce  the  detrimental  effect  of  distortion  Schulmeyer  [43]  conducted  distorted  flow 
experiments  on  a  low  speed  single  stage  compressor  and  showed  that  through  static  restag¬ 
gering  of  inlet  guide  vanes  the  nonuniformity  in  the  velocity  could  be  decreased  by  50%  and 
the  compressor  pressure  rise  could  be  increased  by  5.3%.  However,  the  stalling  mass  flow 
remained  unchanged. 

Epstein  et  al.  [10]  suggested  that  active  control  could  be  used  to  increase  the  stable 
operating  range  of  compression  systems.  Shortly  after  this  researchers  gave  experimental 
evidence  that  proved  the  viability  of  the  idea  by  applying  active  control  to  surge.  We 
refer  to  the  survey  paper  by  Greitzer  et  al.  [18]  for  a  discussion  of  various  experiments  and 
concentrate  here  only  on  the  control  of  rotating  stall.  Working  on  a  low  speed  multi-stage 
compressor  Day  [6]  showed  that  control  scheduling  could  be  used  to  increase  the  stable 
operating  range.  By  sensing  velocity  perturbations  with  an  array  of  hot-wires,  jets  were 
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turned  on  to  blow  in  the  tip  region  for  a  predetermined  time  once  a  stall  cell  was  detected, 
lowering  the  stall  point  by  6%. 

The  first  active  control  of  rotating  stall  was  done  by  Paduano  [37]  and  Paduano  et 
al.  [38]  on  a  low  speed  single  stage  axial  compressor.  By  taking  into  account  that  different 
modes  of  a  compression  system  go  unstable  independently  and  one  after  the  other  as  mass 
flow  decreases,  controllers  were  found,  one  at  a  time,  to  stabilize  the  modes  as  they  went 
unstable.  With  this  approach  the  stable  operating  range  was  increased  by  21%.  The  same 
approach  was  used  by  Haynes  [20]  and  Haynes  et  al.  [21]  on  a  low-speed  three-stage  axial 
compressor  and  increased  the  stable  operating  range  by  7.8%.  In  addition,  Haynes  [21] 
showed  that  unsteady  viscous  effects  could  be  modelled  with  a  simple  first  order  lag  model 
and  identified  the  associated  time  constant  from  experimental  data.  In  both  cases  flow 
perturbations  were  sensed  with  a  set  of  hot-wires  and  movable  inlet  guide  vanes  were  used 
for  actuation. 

By  using  aeromechanics!  feedback  Gysling  [19]  showed  that  is  was  possible  to  lower  the 
stalling  flow  coefficient  of  a  single  stage  low  speed  axial  compressor  by  10%.  Gysling  also 
observed  acoustic  modes  in  the  compression  system  and  showed  that  these  modes  can  lead 
to  instability. 

All  the  control  experiments  discussed  above  were  done  for  circumferentially  uniform 
flow.  The  work  presented  here  extends  active  control  of  rotating  stall  for  the  first  time  to 
the  nonuniform  flow  case.  In  addition,  the  dynamics  of  small  perturbations  in  the  presence 
of  large  distortions  have  never  been  measured  before  and  the  predictive  capabilities  of  the 
Hynes- Greitzer  model  have  not  been  established.  Thus,  the  objective  of  this  research  is  to 
show  that  active  control  can  be  used  to  increase  the  stable  operating  range  in  the  presence 
of  large  distortions.  In  addition,  we  would  like  to  determine  the  predictive  capability  of  the 
model  and  analyze  the  effects  of  distortion  on  the  stability  of  an  experimental  compression 
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system.  We  would  like  to  give  answers  to  the  following  research  questions.  Does  the 
linearized  model  predict  the  dynamics  of  small  perturbations?  Is  the  model  accurate  enough 
for  the  design  of  modern  controllers?  Are  controllers  designed  for  circumferentially  uniform 
flow  effective  in  the  presence  of  distortion?  If  not,  cam  we  find  other  controllers  that  will 
stabilize  the  system? 

The  thesis  is  structured  as  follows.  In  Chapter  2  we  discuss  the  MIT  three-stage  com¬ 
pressor  and  experimental  procedures.  The  modelling  assumptions  are  stated  in  Chapter  3, 
and  the  basic  Hynes-Greitzer  model  is  presented  in  a  state-space  description.  The  model 
is  then  extended  to  include  the  effects  of  unsteady  viscous  effects.  Open  loop  experimental 
results  are  presented  in  Chapter  4.  Active  control  of  rotating  stall  is  discussed  in  Chapter  5. 
A  summary  and  suggestions  for  future  research  are  given  in  Chapter  6. 
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2  Experimental  Compression  System 

All  the  experiments  were  conducted  on  a  low  speed  three-stage  compressor.  In  this  chap¬ 
ter  we  describe  the  compressor,  distortion  generator,  instrumentation,  and  experimental 
procedures. 


2.1  The  MIT  Three-Stage  Compressor 

The  experiments  were  conducted  on  a  low-speed  three-stage  research  compressor  originally 
designed  by  Pratt  and  Whitney.  The  compressor  has  been  used  in  the  past  for  various 
studies  by  Gamier  who  studied  the  behavior  of  small  perturbations  during  stall  inception, 
and  Haynes,  who  modified  the  compressor  for  active  control  by  including  actuator  guide 
vanes  (AGVs)  directly  upstream  of  the  compressor.  Haynes  gives  a  detailed  description  of 
the  design  of  the  AGVs  and  instrumentation. 

To  study  the  effects  of  inlet  total  pressure  distortion  the  rig  was  modified  by  increasing 
the  inlet  duct  length  so  that  a  distortion  generating  screen  could  be  installed  fax  enough 
upstream  of  the  compressor  to  decouple  the  two  components.  Flow  field  perturbations 
induced  by  the  compressor  decay  exponentially  as  e“nlxl  The  upstream  ducts  were  thus 
lengthened  to  2.5  mean  radii  to  decrease  the  interaction  between  the  distortion  screen  and 
compressor.  This  suggestion  is  due  to  Longley  [27].  Separating  the  screen  and  compressor 
by  a  large  distance  is  not  a  limitation  of  the  theory  but  did  simplify  the  experimental 
procedure  because  the  stagnation  pressure  distortion,  which  is  the  input  to  the  flow  field 
computation,  does  not  depend  on  compressor  behavior.  Section  2.2. 

To  facilitate  steady  state  measurements  around  the  annulus  the  distortion  screen  was 
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Table  2.1:  MIT  Three-Stage  Axial  Compressor  Design  Parameters 


Tip  Diameter  (mm.) 

610.0 

Hub-to-Tip  Ratio 

0.88 

Design  Average  Reaction 

0.75 

Design  Flow  Coefficient 

0.59 

Pressure  Rise  Coefficient  (design) 

2.03 

Efficiency  (design) 

84.3% 

Stalling  Flow  Coefficient 

0.460 

mounted  on  a  rotating  ring  driven  by  a  stepper  motor.  This  allows  generation  of  measure¬ 
ments  around  the  annulus  at  any  number  of  points  using  fixed  instrumentation  by  rotating 
the  screen  to  the  desired  location.  The  maximum  speed  at  which  the  screen  can  be  ro¬ 
tated  is  20  rpm.  A  mounting  arrangement  was  also  provided  to  include  a  second  stationary 
screen  directly  upstream  of  the  first  screen  (see  Figure  2.2).  With  the  rotating  and  sta¬ 
tionary  screen  it  is  possible  to  create  distortions  of  varying  length.  A  schematic  layout  of 
the  complete  system  is  shown  in  Figure  2.1,  and  Figure  2.2  shows  a  part  of  the  upstream 
duct  and  distortion  generator.  Tables  2.1  and  2.2  list  the  compressor  design  and  geometric 
parameters.  See  also  Table  3.2  on  page  66. 
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Table  2.2:  MIT  Three-Stage  Axial  Compressor  Parameters. 
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Figure  2.1:  Schematic  of  Three  Stage  Compressor. 
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2.2  Distortion  Generator 

Two  factors  determined  the  choice  of  the  distortion  screens  used  in  the  experiments.  First, 
the  magnitude  of  the  distortion  had  to  be  large  enough  to  introduce  strong  coupling  between 
different  spatial  harmonics  of  perturbations,  so  that  the  dynamic  behavior  with  distortion 
differed  significantly  from  the  undistorted  case.  The  basic  Hynes- Greitzer  model  (discussed 
in  Section  3.3)  was  used  to  determine  the  transfer  functions  between  the  different  harmonics, 
and  distortion  magnitudes  were  chosen  so  that  the  magnitudes  of  the  coupling  transfer 
functions  between  the  first  and  second  harmonics  were  the  same  order  as  that  of  the  first 
harmonic  transfer  function.  (In  the  absence  of  distortion  the  coupling  between  the  first  and 
second  harmonic  is  zero.) 

Modal  analysis  of  the  linearized  dynamics  showed  that  a  distortion  extent  (the  part  of 
the  annulus  that  is  blocked  by  the  screen,  see  Figure  2.3)  of  120  -  240°  has  strong  first, 
second  and  third  harmonics,  satisfying  the  coupling  requirement.  It  will  be  shown  that  an 
extent  of  approximately  120°  is  a  worst  case  distortion  (defined  in  Section  3.7).  Examples  of 
the  harmonic  content  is  given  in  Section  3.6.  Distortion  magnitude  and  extent  are  defined 
in  Figure  2.3. 

Second,  the  distortion  should  be  large  enough  to  give  measurable  changes  in  the  stalling 
flow  coefficient  and  peak  pressure  rise.  The  measured  change  in  stalling  flow  coefficient 
turned  out  to  be  small  for  distortions  up  to  about  one  dynamic  head,  (defined  below)  even 
though  there  is  strong  coupling  present  between  the  harmonics.  It  was  thus  decided  to  carry 
out  experiments  at  two  different  distortion  magnitudes;  one  of  roughly  0.8  dynamic  head, 


Figure  2.3:  Distortion  magnitude  and  extent. 


where  cx  is  the  mean  axial  velocity,  A pt  is  the  drop  in  total  pressure  across  the  screen,  and 
p  the  density  of  air.  The  second  distortion  screen  had  a  magnitude  of  1.9  dynamic  head.  In 
both  cases  the  distortion  extent  was  120°.  Square-shaped  distortions  were  used  as  this  is  a 
standard  used  by  engine  manufacturers.  The  screens  were  designed  using  the  formulations 
given  by  Bruce  [3]  and  Koo  and  James  [24]. 


In  the  following  pressures  will  be  nondimensionalized  by  pu^,  where  itw  is  the  mean 
wheel  speed,  see  Section  3.3.1.  The  distortion  magnitude  defined  above  has  been  nondi¬ 
mensionalized  by  the  inlet  dynamic  head.  The  relation  between  distortion  magnitudes  is 


_  2 d 
dc*  ~  <fp 


(2.2) 


where  d  =  — ?r  is  the  distortion  magnitude  in  dynamic  head  based  on  wheel  speed,  and 

PK 

cp  =  cx/uw  is  the  flow  coefficient. 
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A  parameter  that  is  often  used  to  assess  the  severity  of  inlet  distortion  is  the  DC  (60) 
descriptor  (see  Williams  [49])  defined  by 

DC(60)  =  ^360°  ~  Pt|worst60°  fr,  o\ 

In  an  idealized  case  where  the  static  pressure  is  uniform,  DC(60)=1  corresponds  to  zero- 
velocity  flow  in  a  60°  .  sector  of  the  annulus.  This  is  considered  very  poor  inlet  aero¬ 
dynamics  [49].  Experiments  by  Aulehla  and  Schmitz  [2]  on  the  Tornado  showed  that 
0.15  <  DC(60)  <  0.55  as  the  angle  of  attack  is  varied  over  the  range  3°  to  35°  degrees. 

At  the  respective  stalling  flow  coefficients,  the  0.8  and  1.9  dynamic  head  distortions 
correspond  to  DC(60)=0.53  and  DC(60)=1.31  respectively.  Thus,  the  intensities  of  the  two 
distortions  used  in  the  experiments  represent  typical  and  very  poor  inlet  conditions. 

2.3  Instrumentation 

The  layout  of  most  of  the  instrumentation  is  shewn  in  Figure  2.4.  Steady  state  static 
pressures  were  measured  with  8  upstream  and  8  downstream  wall  taps  on  the  outer  casing 
of  the  compressor.  Total-static  Pitot  probes  were  also  mounted  at  8  equally  spaced  locations 
around  the  annulus  upstream  of  the  AGVs,  at  midspan  these  were  used  to  determine  the 
compressor  pressure  rise  and  calibrate  the  hot-wires.  The  total-static  Pitot  probes  were 
mounted  5°  degrees  away  circumferentially  from  8  of  the  hot-wires  to  ensure  that  they  did 
not  interfere  with  the  flow  at  the  hot-wire  tips. 

Two  total-static  Pitot  probes  were  mounted  downstream  of  the  distortion  screen  (not 
visible  in  Figure  2.4).  All  the  pressure  probes  were  connected  to  a  48-channel  Scanivalve 
that  was  controlled  by  a  VAXstation. 
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Axial  Station 


41 


Figure  2.4:  Sensor  Locations. 


Mass  flow  was  measured  by  an  orifice  plate  fax  downstream  of  the  throttle  (see  Fig¬ 
ure  2.1).  The  pressure  drop  across  the  orifice  plate  was  measured  by  an  MKS  Baratron.  The 
standard  deviation  of  the  flow  coefficient  measurements  reported  here  is  less  than  0.08%. 

Rotor  rotation  speed  was  measured  with  a  magnetic  pick-up  and  60  tooth  gear  mounted 
on  the  drive  shaft. 

Flow  velocities  were  measured  at  midspan  upstream  of  the  AG  Vs  with  16  equally  spaced 
Dantec  55-P11  hot-wires  connected  to  Dantec  56C17  bridges  and  CTA  56C01  units.  The 
hot-wires  were  calibrated  before  each  experiment  over  the  range  6  m/s  to  60  m / s  (equivalent 
flow  coefficient:  0.08  -  0.8);  the  calibration  procedure  is  discussed  in  Section  2.4.1.  The  hot¬ 
wire  signals  were  filtered  with  fourth  order  Bessel  filters  before  they  were  sampled  at  500  Hz. 
The  filter  cut-off  frequency  was  set  to  200  Hz.  This  cut-off  frequency  was  low  enough  to  give 
approximately  80  dB  attenuation  at  the  rotor  blade  passing  frequency.  Analogic  HS  DAS-16 
16-bit  A/D  converters  with  AMUX-64-X  64  channel  multiplexer  were  used  to  discretize  the 
hot-wire  signals.  The  resolution  of  the  A/D  was  better  than  0.01  m/s. 

The  12  AGV  servo  motors  were  controlled  by  DMC  430  servo  motion  controllers.  These 
controllers  ran  at  2  kHz  sampling  rate  asynchronously  from  the  main  control  loop.  AGV 
deflections  were  measured  by  4096  pulse  per  revolution  incremental  shaft  encoders. 


2.4  Experimental  Procedures 

In  this  section  we  outline  the  hot-wire  calibration  and  measurement  of  the  stalling  flow 
coefficient  procedures. 
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2.4.1  Hot-wire  Calibration 


A  new  hot-wire  calibration  procedure  was  developed  to  deal  with  the  circumferentially 
nonuniform  flow  around  the  annulus.  Eight  total-static  Pitot  probes  were  used  to  measure 
the  velocities  around  the  annulus.  There  are  large  gradients  in  the  steady  flow  profile  at 
the  edges  of  the  distorted  region  and  it  was  necessary  to  include  the  5°  separation  between 
the  hot-wires  and  total-static  probes  in  the  calibration  procedure. 

The  calibration  procedure  was  basically  to  measure  the  flow  velocity  with  the  total- 
static  probes  and  then  rotate  the  screen  so  that  the  measured  velocities  are  aligned  with 
one  set  of  eight  hot-wires.  The  actual  procedure  used  is  more  involved  and  is  explained 
with  the  aid  of  Figure  2.5  and  the  steps  below.  In  the  figure  we  assume  the  first  hot-wire  is 
at  0°  and  the  first  total-static  probe  at  5°.  Recall  that  the  hot-wires  are  separated  by  22.5° 
so  the  angle  between  a  total-static  probe  and  the  next  hot-wire  is  17.5°. 

1.  Select  a  mass  flow  and  execute  the  following  four  steps. 

(a)  Measure  the  velocities  with  the  8  total-static  probes.  The  measured  velocities 
are  indicated  by  circles  in  the  graph  at  the  top  of  Figure  2.5. 

(b)  Rotate  the  screen  through  5°  and  measure  all  the  hot-wire  voltages.  The  voltages 
of  the  hot-wires  at  0,45°,...  correspond  to  the  velocities  obtained  from  the 
previous  step.  This  is  shown  in  the  second  graph  from  the  top  in  the  figure. 

(c)  Rotate  the  screen  through  17.5°  and  measure  the  velocities  with  the  total-static 
probes.  These  velocities  correspond  to  the  hot-wire  voltages  at  22.5, 67.5°, .. .  in 
the  previous  step.  These  velocities  are  indicated  by  stars  in  the  third  graph  of 
the  figure. 

(d)  Rotate  the  screen  through  5°  and  measure  all  the  hot-wire  voltages.  The  voltages 
of  the  hot-wires  at  0,45°,...  correspond  to  the  velocities  obtained  from  the 
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previous  step  —  these  velocities  axe  indicated  by  stars  in  the  bottom  graph  in 
the  figure.  In  addition,  the  velocities  measured  in  step  1(a)  are  now  aligned  with 
the  hot-wires  at  22.5, 67.5°, . . .  and  is  indicated  by  the  circles  in  the  graph  at  the 
bottom  of  the  figure. 

Steps  l(a)-l(d)  give  two  calibration  points  for  each  hot-wire. 

2.  Steps  l(a)-l(d)  are  repeated  three  more  times  at  90°  intervals,  giving  a  total  of  eight 
calibration  points  for  each  hot-wire.  This  step  uses  the  fact  that  the  flow  is  nonuniform 
axound  the  annulus  and,  instead  of  changing  the  mass  flow,  we  move  to  a  different 
point  on  the  profile  to  get  a  different  velocity. 

3.  Steps  1  and  2  above  are  repeated  at  two  more  flow  coefficients  to  give  a  total  of 
24  velocities  for  each  hot-wire.  The  three  different  flow  coefficients  axe  chosen  so 
that  the  absolute  minimum  to  absolute  maximum  velocity  ranges  from  approximately 
6  m/s  to  60  m/s.  This  range  covered  all  the  velocities  that  were  encountered  during 
experiments.  King’s  law  [14]  was  used  to  determine  a  relation  between  the  hot-wire 
voltages  and  corresponding  velocities. 

2.4.2  Stalling  Flow  Coefficient  Measurement 

It  was  found  that  the  stall  point  was  sensitive  to  the  rate  at  which  the  throttle  was  closed. 
It  was  thus  necessary  to  develop  a  procedure  by  which  the  stall  point  could  be  determined 
in  a  consistent  way. 

With  a  rough  estimate  of  the  stalling  flow  coefficient  available,  the  throttle  was  closed 
slowly  until  the  mean  flow  was  about  3%  above  the  stall  point.  From  this  point  on  the  throt¬ 
tle  was  closed  at  the  minimum  possible  rate  until  the  expected  stall  point  was  reached  — 
this  typically  took  about  two  minutes.  The  throttle  was  then  stopped.  The  system  had  to 
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remain  stable  for  at  least  two  minutes  at  this  flow  coefficient  before  it  was  accepted  as  the 
possible  minimum  point.  If  the  compressor  was  still  stable  after  two  minutes,  the  throttle 
was  again  closed  at  minimum  rate  until  the  flow  coefficient  was  reduced  by  approximately 
0.2%.  Again  the  system  had  to  remain  stable  for  at  least  two  minutes  at  this  flow  coefficient. 
This  last  step  was  repeated  until  a  stall  point  was  reached.  After  a  stall  the  throttle  was 
opened  until  the  compressor  returned  to  unstalled  operation  and  then  the  procedure  was 
repeated  from  the  beginning.  If  the  same  flow  coefficient  was  reached  three  times  it  was 
accepted  as  the  stall  point  for  the  particular  distortion/controller  configuration.  In  most 
cases  experiments  were  repeated  on  different  days  and  the  repeatability  of  the  stalling  flow 
coefficients  was  better  than  0.2%. 
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3  Compression  System  Model 


In  this  chapter  we  develop  a  complete  model  for  a  compressor  operating  in  the  presence  of  an 
inlet  circumferential  total  pressure  distortion.  The  goal  is  to  obtain  a  linearized  state  space 
description  of  the  behavior  of  small  perturbations  about  a  circumferentially  nonuniform 
steady  flow.  We  start  by  stating  the  ideas  that  underlie  the  model.  Most  variables  are 
periodic  and  it  is  useful  to  express  them  as  Fourier  series.  Manipulating  Fourier  series 
coefficients  is  convenient  if  we  vectorize  the  coefficients;  the  mathematical  tools  for  doing 
so  are  given  in  Section  3.2.  Derivation  of  the  linearized  dynamics  is  done  in  two  steps. 
In  the  first  we  determine  a  steady  operating  point  (mean  mass  flow  or  throttle  setting) 
by  solving  one  or  more  sets  of  nonlinear  paxticil  differential  equations.  We  then  linearize 
the  nonlinear  partial  differential  equations  about  the  operating  point.  These  processes  are 
described,  as  is  the  extension  of  the  basic  analysis  to  include  a  simple  model  of  unsteady 
blade  row  response.  The  model  is  then  used  to  compare  uniform  and  nonuniform  flow  mode 
shapes,  and  explain  the  control  theoretic  difference  between  uniform  and  distorted  flow. 


3.1  Assumptions 

In  this  section  we  summarize  the  relevant  modelling  assumptions.  Detailed  discussion  is 
given  by  Haynes  [20]. 

The  flow  is  incompressible.  In  the  experiment  mean  blade  Mach  number  was  0.2  and 
the  modal  wave  Mach  number  0.1. 

The  flow  is  two  dimensional  so  that  radially  averaged  description  is  used.  Haynes 
showed  that  this  was  a  good  approximation  for  high  hub-to-tip  ratio  compressors.  The 
hub-to-tip  ratio  of  the  compressor  used  is  0.88. 
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Viscous  effects  axe  negligible  outside  the  blade  rows.  This  assumption  is  valid  if  inertial 
effects  are  much  larger  than  viscous  effects.  In  the  experiment  the  Reynolds  number  for  the 
nth  harmonic  was  on  the  order  of  6  •  105/n  so  that  this  assumption  is  justified. 

The  compressor  can  be  modelled  as  a  semi-actuator  disk.  The  compressor  is  the  most 
complex  component  in  the  system.  This  assumption  allows  us  to  describe  the  pressure 
rise  across  the  compressor  in  a  simple  way  without  describing  the  detail  of  the  flow  in  the 
individual  blade  rows.  The  details  of  the  assumption  will  be  discussed  when  we  present 
expressions  for  the  pressure  rise  across  the  compressor. 

It  is  important  to  note  that  two  assumptions  used  in  previous  models  for  the  control 
of  rotating  stall  have  not  been  included.  The  first  is  the  assumption  that  the  steady  flow  is 
circumferentially  uniform.  We  are  specifically  interested  here  in  the  behavior  of  the  com¬ 
pressor  in  the  presence  of  circumferentially  nonuniform  steady  flow.  The  second  assumption 
that  has  been  dropped  is  the  unimportance  of  overall  system  (surge)  dynamics.  When  the 
steady  flow  is  circumferentially  uniform  and  the  velocity  perturbations  are  small  there  is 
no  coupling  between  the  surge  and  rotating  stall  modes  so  that  the  surge  dynamics  can 
be  ignored.  When  the  steady  flow  is  circumferentially  nonuniform,  however,  strong  cou¬ 
pling  can  exist  between  the  surge-like  and  rotating  stall  modes.  It  is  necessary  to  include  a 
description  of  the  overall  compression  system  as  well. 


3.2  Spatial  Fourier  Series 

Most  variables  in  the  model  will  be  functions  of  0  and  it  is  convenient  to  express  them  as 
Fourier  series,  that  is,  any  function  f(9,  r)  will  be  written  as 

/(0,t)  =  fo(T)  +  ^r(fcn(T)cos(ne)  +  fan(T)sin{nO)).  (3.1) 

n>0 
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The  function  /0(r)  will  be  referred  to  as  the  zeroth  harmonic  and  will  always  be  included 
in  the  expansion.  Harmonics  one,  two,  . . .  will  be  referred  to  as  higher  harmonics.  Note 
that  we  allow  for  time  varying  Fourier  coefficients. 

The  development  of  the  model  is  simplified  considerably  if  we  use  vectors  to  represent 
Fourier  series  coefficient's.  To  do  so,  we  stack  the  Fourier  coefficients  in  a  vector  f  as  follows 


/  —  [/o>  /cij/sii  /c2>/s2)  •  •  ■  > /cn> /sn] 


Multiplication  of  Fourier  series  with  known  coefficients  by  Fourier  series  with  unknown 
coefficients  occurs  frequently.  Let  g{8)  be  a  known  Fourier  series;  we  would  like  to  write 
the  coefficient  vector  p  of  the  product  p{8)  =  g{8)f{8)  in  terms  of  the  coefficient  vector 
/.  Multiplication  of  two  functions  in  the  spatial  domain  corresponds  to  convolution  in 
the  Fourier  domain,  thus  the  coefficients  of  the  product  can  be  computed  by  a  matrix 
multiplication.  Let  F(g)  be  the  Fourier  convolution  matrix  corresponding  to  g(8),  then  we 
can  write 

p(8)=g(8)f(8)  >  p  =  F(g)f.  (3.3) 

The  construction  of  F(g)  follows  from  standard  trigonometric  identities  —  the  detail  is 
presented  in  Appendix  A.  We  will  assume  that  the  dimensions  of  the  Fourier  convolution 
matrices  have  been  selected  appropriately  so  that  all  multiplications  and  additions  are  well 
defined  in  the  equations  that  follow. 

Derivatives  with  respect  to  8  can  also  be  written  as  a  matrix  multiplication.  With  f{8) 
defined  as  before,  let 

m  =  %  0.4) 

=  0  +  ^(n/sn(r)  cos(n0)  -  nfcn{r)  sin(n0)).  (3.5) 

n>0 
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Vectorizing  this  we  get 


9  =  Dgf 


(3-6) 


where 


Ds  =  diag 


- 1 

o 

0  1 

> 

0  2 

)  *  1  1 

o 

r—\ 

1 

_ 1 

1 

1 

to 

L_°_ 

(3.7) 


Equations  (3.3)  and  (3.6)  allow  us  to  vectorize  all  the  expressions  in  the  model. 


3.3  Model  With  Quasi-Steady  Viscous  Flow  Effects 

We  will  study  the  compression  system  shown  schematically  in  Figure  3.1.  The  system  con¬ 
sists  of  a  distortion  screen,  long  upstream  duct,  actuator  guide  vanes  (AGVs),  a  compressor, 
downstream  duct,  plenum,  and  throttle.  The  experimental  compressor  has  a  row  of  inlet 
guide  vanes  directly  upstream  of  the  AGVs  (see  Figure  2.4)  that  is  not  shown  in  this  figure 
but  will  be  accounted  for  in  the  modelling.  Also  shown  in  the  figure  is  the  axis  system  that 
will  be  used.  The  origin  of  the  axial  variable  is  taken  at  the  face  of  the  compressor. 

3.3.1  Introduction 

The  following  parameters  will  be  used  to  nondimensionalize  the  various  quantities. 

f  =  mean  rotor  radius 

uw  =  mean  wheel  speed 

pv%,  =  dynamic  head  based  on  wheel  speed, 

p  —  density  of  air  under  operating  conditions. 
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Distortion 


AG  Vs  Compressor 


Plenum 


Figure  3.1:  Compression  System. 

The  following  nondimensional  independent  variaDies  axe  used. 

9  —  angle  around  annulus,  positive  in  direction  of  rotor  rotation. 

r  —  nondimensional  time  =  (time)uw/f 

x  =  nondimensional  axial  position,  origin  is  at  compressor  face  (see  Figure  3.1) 

The  state  of  the  system  is  described  in  terms  of  the  axial  and  circumferential  flow 
velocities,  and  total  and  static  pressure.  The  model  is  presented  in  nondimensional  form 
only  and,  to  simplify  notation,  the  nondimensionalization  will  not  be  shown  explicitly,  e.g., 
we  will  write  p  instead  of  p/{pu ^).  The  main  variables  are 
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Ps(0,T,x) 

Pt{0,T,x) 

4>{9,t,x) 

7{9,t) 


(static  pressure  -  patmosPhere)/(/”4) 

(total  pressure  -  patmosphere)/(/™w) 

flow  coefficient  =  (axial  velocity) /uw 

AGV  deflection  angle;  we  assume  a  continuum  of  blades. 


The  dependence  on  the  independent  variables  will  often  be  omitted. 

The  semi-actuator  disk  model  uses  the  geometry  of  the  compressor  and  the  axisym- 
metric  compressor  characteristic.  The  parameters  are 


ip(4>,'r)  = 


M4>)  = 


Pr 


Ps 


Pi 


Pc 


Pa. 


X 


m 


nondimensional  steady  state  axisymmetric  total-to- 
static  pressure  rise  across  the  compressor, 
nondimensional  isentropic  total-to-static  pressure  rise 
which  follows  from  the  torque  work  done  on  the  air. 
rotor  fluid  inertia  =  ^  c~sjr;j  w^ere  j  is  the  nondi¬ 
mensional  rotor  blade  chord  at  midspan  for  the  jth 
stage,  and  7 Tj  is  the  corresponding  rotor  stagger  angle, 
same  as  nT  but  for  stators. 

same  as  (j.r  but  for  inlet  guide  vanes.  The  inlet  guide 
vanes  are  not  shown  in  Figure  3.1,  see  Figure  2.4  in¬ 
stead. 

pT  +  ps  +  Pi  =  compressor  fluid  inertia. 

AGV  fluid  inertia  =  cost*)  w^ere  Tso  is  the  circumfer¬ 
ential  mean  AGV  deflection. 

nondimensional  upstream  axial  location  where  all  ve¬ 
locity  and  pressures  are  measured.  Note  that  xm  <  0. 
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Steady  state  values  will  be  indicated  by  a  San  Serif  subscript  s  to  distinguish  them 
from  the  Roman  subscript  s  that  denotes  stators.  Station  numbers  as  defined  in  Figure  3.1 
will  be  denoted  by  parenthesized  superscripts.  All  variables  will  be  written  as  the  sum  of  a 
steady  state  value  plus  a  time  varying  perturbation;  for  example,  the  flow  at  the  compressor 
face  is 

<t>{3)(9,T)=<t>™(6)  +  6<f>w(d,T)  (3.8) 

where  5  indicates  perturbation. 

3.3.2  Pressure  Balances 

To  derive  the  equations  that  express  the  flow  in  the  compression  system,  pressure  balances 
will  be  written  for  each  component  and  then  combined  in  an  overall  pressure  balance  for 
the  complete  system. 

AGV  pressure  balance 

The  flow  coefficient  and  pressure  balance  relations  for  the  AGV  operating  under  nonuniform 
steady  flow  are  derived  in  Appendix  B  and  are  given  by 

8<t>m  =  (1  +  +  M*^7  (3-9) 

=  /*ty(3)+M 7  (3-10) 

<Xp[3)  -  =  -Mad  +  y||w(3)  -  y&|*7  (3.11) 

=  P^4>{3)  +Pi5 7- 


(3.12) 


The  expressions 


U  =  1  +  ^ 

(3.13) 

tn 

-e- 

II 

< 

(3.14) 

defined  in  Equation  (3.10)  are  independent  of  time  and  depend  on 

the  steady  state  AGV 

deflection  qs  and  flow  <f>s  respectively  so  that  they  are  functions  of  9.  The  same  holds  true 

for 

p^  =  -^(l  +  T—) 

(3.15) 

n  d 

Pi~  2*%de 

(3.16) 

defined  in  Equation  (3.12).  Vectorizing  Equations  (3.10)  and  (3.12)  using  Equations  (3.3) 

and  (3.6)  we  get 

54>i2)  =  F ^  +  F^S 7 

(3.17) 

6(p[i)-p[7))  =  P^W  + 

(3.18) 

where 

(3.19) 

1 Ry  ==  ^(/^a^s 

(3.20) 

(3.21) 

•P-r  =  -§F(<f>s)Dg. 

(3.22) 

Upstream  flow  field 

By  taking  the  first  integral  of  the  axial  momentum  equation  and  nondimensionalizing  we 
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obtain  [22] 


£(Pt2)  -  Pt])  =  4^’  -  “  (^cn  cos(n^)  +  <5^  sin(n0))  .  (3.23) 

1  n> 0  71 

The  equation  is  linear  in  the  coefficients  of  5ft2)  so  we  can  vectorize  this  by  constructing  a 
matrix 


A21  =  —  diag 


I  ii  II 

IV  '  '  2’  2’ 


1  1 

3’  3  ’  ’ 


(3.24) 


so  that 


6(j>?)-p?))  =  An6<t> 


~(2) 


(3.25). 


We  will  work  in  terms  of  6<pm,  the  flow  through  the  compressor.  Substituting  the  time 
derivative  of  Equation  (3.17)  into  Equation  (3.25)  gives 

~(3) 

6(p™  -  p'l))  =  +  A21F76 7.  (3-26) 


Compressor  Pressure  Balance 


The  total-to-static  pressure  rise  across  the  compressor  has  been  developed  by  Hynes  and 
Greitzer  [22] 

Ps4)  ~  Pt3>  =  W,  7)  “  ~  P<^<3)  (3-27) 

where  ^(<£,7)  IS  the  steady  inlet-total  to  exit-static  pressure  rise  in  axisymmetric  flow,  /xr 
is  the  rotor  fluid  inertia,  and  Me  =  Mr  +  Ms  +  Mi  IS  the  total  fluid  inertia  in  the  compressor. 
Linearizing  this  equation  about  the  nonuniform  steady  flow  4>s  we  get 


«(rf>  -?!”)  = 


dtp  d 
d$~^rde 


6<f>(3)  +  -  Pc^(3)- 

07 


(3.28) 
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Recall  that  ^  and  ^  are  functions  of  0  Vectorising  the  last  equation  gives 

4®*’  -pf)  - 


6<t>{3)  -  fxc6(f>{  +  F{^-)8^. 


(3.29) 


Downstream  Ducts 

Here  we  assume  that  the  steady  state  compressor  exit  flow  through  a  row  of  high  solidity 
stators  is  axial  and  that  there  are  no  downstream  obstructions.  Under  these  conditions  the 
static  pressure  is  circumferentially  uniform  and  will  be  the  same  as  the  plenum  pressure  [22] 

P?’-P ?’•  (3-30) 


Using  mass  conservation  through  the  compressor  and  the  high  solidity  exit  condition,  Hynes 
and  Greitzer  [22]  have  shown  that  the  downstream  static  pressure  perturbations  satisfy 


-  P(s4))  =  -j-4?  -  E  £  (*$8  co<nd)  +  4%  sin(n0))  •  (3-3!) 


i 


n>0 


We  vectorize  this  equation  in  the  same  way  as  Equation  (3.23)  to  get 


6(p(s5) 


p™)  =  ASM 


.(3) 


(3.32) 


where 


A54  =  —  diag 


I  ii  II  H 

14  l  l  o  o 


(3.33) 


Plenum  Dynamics 

The  nonuniform  steady  flow  introduces  coupling  between  the  different  harmonics,  including 
the  zeroth  harmonic,  and  it  is  thus  necessary  to  include  pressure  perturbations  in  the  plenum 
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in  the  description.  Assuming  that  the  pressure  and  density  changes  in  the  plenum  can  be 
related  isentropically  and  that  the  flow  through  the  throttle  can  be  described  by 

(3-34) 

where  kt  is  the  throttle  constant,  we  obtain 


cfrso 


1 

462/t 


6<t>? 


(3.35) 


where 


Zt  =  +  lA  +  +  fxc  (effective  compressor  length)  (3.36) 

6  _«w  /Plenum  volume  ^ 

G  2 us  V  Duct  area  x  lt 

and  ua  is  the  speed  of  sound.  Note  that  we  assume  Sp(s6)  to  be  nonzero  and  thus  it  acts  as 
an  external  forcing  function.  Only  the  last  term  on  the  right  hand  side  of  Equation  (3.35) 
needs  vectorization.  All  we  need  to  do  here  is  construct  a  matrix  that  selects  the  zeroth 
harmonic  of  <5^(3).  The  matrix  has  a  single  row  with  the  first  element  equal  to  one  and  the 
rest  of  the  elements  zero.  Thus 


=  *-  M"  -  M”)  +  -W” 

4&QZt^t0so 
As+  =  ^462/?  °’ 


(3.38) 

(3.39) 


This  completes  the  pressure  balances  for  the  compression  system. 

The  overall  pressure  balance  for  the  perturbations  is  given  by 

6(p ?  -  =  (M6)  -  6p?)  +  (M5)  ~  ¥s4))  +  (H4)  ~  ¥*) 

+  (SKm-SKm)  +  (S*m-rtl))-  (3-4°) 
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Substituting  from  the  corresponding  equations  and  simplifying  we  get 


'  ^(3)' 

8(f) 

a-a,5 

64>i3) 

M<t> 

-l/(46&Jt*t&o)_ 

6p?\ 

+ 


A^A^+Py] 


0 


+ 


0 

Sp[» 

1/(4  b2altkt^s0)_ 

=  A8x  +  PJ7  +  E87  +  G[5p[l\  SpWf 


(3.41) 


(3.42) 


where 


A<p  =  Hcl  ~  A21F4  —  —  A54  (3.43) 

A*5  =  [l,  0,  0,  ...]t  (3.44) 

and  I  is  the  identity  matrix.  We  will  denote  the  state  vector  [ 5(j> ,  <fy?s5)]T  by  Sx. 

Note  that  we  assumed  the  total  pressure  perturbation  is  nonzero.  We  thus  have 
two  external  sources,  Sp{1}  and  8pi6\  that  force  the  system.  These  sources,  called  process 
noise  in  the  control  and  estimation  literature,  will  be  used  to  study  the  effect  of  external 
noise  on  the  system. 

~(2m) 

When  we  compare  with  experiments  the  velocity  measurements  8(f)  are  taken  up¬ 
stream  of  the  AGVs.  We  must  therefore  take  the  decaying  potential  field  into  account;  the 
nth  harmonic  decays  with  enim  (x  is  negative  upstream).  The  state  vector  is  defined  in 
terms  of  8<f>(3)  so  we  also  need  the  relation  between  8<f>m  and  8<f>(3)  given  in  Equation  (3.17). 
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Thus 


E^  =  diag[l,  eIm,eIm,  e2l”,e2*m,. . .]  (3.45) 

54>{Xm)  =  E^64>(2)  (3.46) 

=  e(* m) o  +  &Xm)FySi  (3-47) 

=  Ctf*  +  £>£7.  (3.48) 


The  zero  at  the  end  of  the  matrix  in  Equation  (3.47)  multiplies  the  plenum  pressure  per¬ 
turbations,  the  last  state  in  the  state  vector  6x. 

Equation  (3.42)  contains  the  derivative  of  the  control  signal  with  respect  to  time  so 
that  it  is  not  in  the  usual  state  space  form.  It  can  be  converted  to  the  usual  form  with  the 
following  transformation.  Let 

8z  =  5x-  Efn.  (3.49) 

Applying  this  transformation  to  Equations  (3.42)  and  (3.48)we  get 

8z  =  ASz  +  (B  +  AE)8:y  +  G[<5p(x),  <5p'6)]T  (3.50) 

50(lm)  =  CSz  +  (D  +  AE)8i  (3.51) 

that  has  the  desired  form.  We  note  that  the  final  state  space  description  given  by  Equa¬ 
tions  (3.50)  and  (3.51)  is  input-output  equivalent  to  Equations  (3.42)  and  (3.48)  but  not 
state  equivalent.  We  can  recover  the  original  state  vector  8x  from  the  relation  in  Equa¬ 
tion  (3.49)  if  it  is  needed.  However,  for  the  homogeneous  system  8j  =  0  so  that  8z  =  8x 
and  we  can  use  Equation  (3.50)  for  eigenstructure  analysis.  The  model  just  derived  will  be 
referred  to  as  the  basic  Hynes-Greitzer  model. 

Before  we  can  compute  the  matrices  of  the  state  space  model  we  need  to  find  the 
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nonuniform  steady  flow  <f>s.  The  equation  that  needs  to  be  solved  and  solution  method  are 
discussed  in  the  next  section. 

3.3.3  Steady  State 

The  equation  describing  the  steady  state  is  obtained  by  writing  a  pressure  balance  for  the 
system  as  we  did  in  the  previous  section,  but  this  time  the  pressure  balance  is  for  the  steady 
variables.  Because  we  are  interested  in  the  steady  state  solution  all  derivatives  with  respect 
to  time  can  be  set  to  zero  to  give 

P?  ~  P?  =  -*(<M  +  fh^  +  \k A  (3.52) 

The  equation  can  be  vectorized  as  before;  the  detail  is  omitted  here.  Once  we  have  decided 
on  an  operating  point  (by  specifying  either  the  mean  flow  <j)s0  or  throttle  constant  fct),  and 
the  shape,  magnitude,  and  extent  of  the  total  pressure  distortion  p[x)(0),  Equation  (3.52) 
can  be  solved  for  the  remaining  unknown  variables.  Assume,  as  would  be  the  case  in  an 
experiment,  we  have  specified  kt  so  that  all  the  coefficients  of  <ps  are  unknown.  In  general  we 
will  need  an  infinite  number  of  harmonics  for  <f>s.  For  practical  purposes  we  approximate  the 
solution  with  a  finite  number  of  harmonics  so  that  the  right  hand  side  of  Equation  (3.52) 
will  not  equal  the  left  hand  side.  By  defining  the  error  in  the  approximation  as  the  difference 
between  the  left  and  right  hand  side  of  Equation  (3.52)  we  can  find  a  least  squares  solution 
for  (j)s. 

We  are  mainly  interested  in  computing  the  linearized  dynamics  although  steady  state 
quantities  are  also  important.  We  will  use  the  relative  error  in  the  dominant  poles  and  zeros 
as  a  measure  to  determine  the  number  of  harmonics  in  the  steady  velocity  <f>s(0)  and  number 
harmonics  needed  for  the  linearized  dynamics  5</>.  The  “exact”  values  of  the  dominant  poles 
and  zeros  are  defined  as  those  values  obtained  with  a  large  number  of  harmonics.  Once 
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we  have  determined  the  number  of  harmonics  needed  to  obtain  the  “exact”  values,  we 
determine  the  number  of  harmonics  needed  to  give  a  desired  relative  error,  for  example 
0.1%.  By  its  very  nature  the  procedure  to  determine  the  number  of  harmonics  is  iterative. 

The  actual  number  of  harmonics  needed  for  a  good  approximation  is  quite  low.  Through 
simulations  it  was  found  that  if  we  want  to  compute  the  first  n  modes  (that  is,  the  eigenvalues 
(or  poles)  and  corresponding  eigenvectors)  of  the  linearized  dynamics  accurately  we  need 
approximately  8 n  harmonics  for  <f>s.  Once  we  have  solved  for  <ps  we  keep  only  the  first  4 n 
harmonics  and  use  this  truncated  Fourier  series  to  compute  the  linearized  model.  We  need 
to  start  with  about  twice  as  many  harmonics  for  <f>s  as  we  are  interested  in  to  ensure  that 
the  4n  harmonics  that  we  need  are  computed  accurately.  This  can  be  explained  as  follows. 
If  we  multiply  cos p6  by  cos  q9  we  get  products  of  the  form  cos (p  +  q)6  and  cos(p  —  q)9  so 
that  higher  harmonics  effects  the  lower  harmonics. 

Similarly,  we  need  approximately  4 n  harmonics  for  the  perturbations  <5<£(3)  in  the  state 
vector.  We  will  see  later  that  only  the  first  four  modes  are  needed  to  model  the  system  for 
frequencies  up  to  three  rotor  revolutions  so  that  we  need  approximately  16  harmonics  for 
5<t>{3)  and  initially  32  harmonics  (plus  the  zeroth  harmonic)  for  <f>%.  Each  harmonic  needs 
two  coefficients  so  that  Equation  (3.52)  represents  41  coupled  differential  equations  and  the 
order  of  the  state  space  model  is  22.  Solving  the  steady  flow  and  computing  the  state  space 
model  takes  about  2  minutes  on  an  IBM  RS6000  work  station.  The  numbers  quoted  are 
conservative  and  apply  to  a  square  wave  distortion.  For  distortions  without  discontinuities, 
e.g.,  a  cosine,  all  the  computations  can  be  done  with  4n  harmonics. 

All  the  experiments  indicated  that  (for  this  particular  compressor)  four  modes  were 
enough  to  capture  the  dynamics  up  to  three  rotor  revolutions.  It  thus  is  possible  to  reduce 
the  order  of  the  final  state  space  model  from  22  down  to  eight  for  controller  design.  However, 
to  get  to  this  low  number  of  states  we  need  to  start  out  with  a  rather  large  number  of 
harmonics. 
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3.4  A  Simple  Approximation  for  Unsteady  Viscous  Behavior 


In  the  semi-actuator  disk  model  an  instantaneous  change  in  flow  through  the  compressor 
results  in  an  instantaneous  change  in  pressure  rise.  Experiments  by  Nagano  [35]  and  Maz- 
zawy  [32]  have  shown  that  total  pressure  changes  lags  a  change  in  the  flow  and  that  this 
lag  can  be  approximated  quite  well  by  a  simple  first  order  system  with  time  constant  one 
to  two  times  the  flow-through  time.  For  simplicity,  all  the  lags  in  the  rotors  will  be  lumped 
together.  Similarly,  the  lags  in  the  stators  will  also  be  lumped  together.  Haynes  et  al  [21] 
used  this  approach  with  great  success.  We  should  note  here  that  Longley  [25]  has  found  that 
in  some  compressors  it  is  more  important  to  model  deviations  rather  than  losses.  Haynes 
et  al  [21]  found  deviation  lags  not  to  be  important  for  this  compressor.  The  absence  of 
deviations  here  does  not  limit  the  generality  of  the  model  or  the  results  and  can  be  added 
if  the  need  arises. 

The  additional  dynamics  can  be  included  in  the  basic  Hynes-Greitzer  model  by  malt¬ 
ing  two  changes.  First,  we  modify  the  expression  that  describes  the  pressure  rise  across 
the  compressor  to  include  the  effect  of  the  lags.  Then  we  add  the  additional  differential 
equations  to  the  state  space  model,  see  Haynes  et  al  [21]. 

The  assumption  is  that  the  total  pressure  losses  can  be  modelled  by 

l  =  h  +  k  (3.53) 

=  rpi-ip  (3.54) 


where 

l  =  total  pressure  loss  in  rotors  and  stators 
=  /o  +  ^  kn  cos(n^)  +  hn  sin  (n6)  (3.55) 

n>0 
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ipi  is  the  ideal  inlet-total  to  exit-static  pressure  rise,  and  lT  and  la  axe  the  losses  in  rotors  and 
stators  respectively.  It  is  further  assumed  that  the  contribution  of  the  losses  in  the  rotors 
and  stators  is  in  proportion  to  the  reaction  r  of  the  compressor  blading 


lT  =  rl 

(3.56) 

ls  =  (l-r)l. 

(3.57) 

With  these  definitions  we  can  write 

iP  =  ipi-lI-  la. 

(3.58) 

The  second  change  is  the  addition  of  differential  equations  for  rotor  and  stator  losses.  The 
lags  for  the  rotor  and  stator  losses  can  be  modelled  by  first  order  systems 


/dir  ,  dlr\ _  i  ,  i 

Tr^+ae)~~lr  +  ln 
dls  _  ,  ,  , 

T*dr  ~  Is  +  hs- 


(3.59) 

(3.60) 


The  time  constants  rr  and  rs  are  proportional  to  the  flow-through  times  in  the  rotors  and 
stators  respectively 


T’m  =  T,m 

(3.61) 

T’{e)  =  T,m 

(3.62) 

where  bT  and  bs  are  the  mean  rotor  and  stator  lengths  respectively,  and  rf  is  a  constant 
typically  in  the  range  1.0-1. 5.  Linearizing  Equations  (3.59)  -  (3.62)  about  the  steady  state 
values  <f> s,  lIS,  and  las  we  get 


d8lt 

dr 

dSk 

dr 


/  1  .  d  .  1  Sl/S  r,(3) 

+  oZ>5lr  + - 7Z°V 

Tf  86  Tr  5<p 


--a.  + 

Ts  Ts  0(p 


(3.63) 

(3.64) 
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These  equations  can  be  vectorized  and  included  in  the  basic  analysis.  The  results  are 
summarized  in  Table  3.1  on  page  65.  This  model  will  be  referred  to  as  the  extended  Hynes- 
Greitzer  model  (XHG). 

The  parameters  of  the  three-stage  compressor  is  listed  in  Table  3.2  on  page  66. 


Steady  State 

Inclusion  of  the  unsteady  loss  response  also  changes  the  equation  for  the  steady  flow  and  we 
discuss  the  necessary  changes  here.  As  before  we  need  to  modify  the  pressure  rise  expression 
and  add  a  set  of  differential  equations.  The  expression  for  the  pressure  rise  modification 
has  already  been  given  in  Equation  (3.58).  In  steady  state  Equation  (3.60)  simplifies  to  a 
trivial  identity  and  in  Equation  (3.59)  we  simply  put  the  derivative  with  respect  to  time 
equal  to  zero.  The  final  set  of  differential  equations  that  must  be  solved  are  also  shown  in 
Table  3.1. 

To  solve  the  steady  flow  equations  we  need  initial  estimates  for  <ps  and  ls.  A  good  initial 
guess  is  the  steady  state  solution  from  the  basic  model,  with  zero  as  an  initial  estimate  for 
the  losses.  Furthermore,  it  was  found  that  the  number  of  harmonics  needed  in  all  the 
calculations  can  be  reduced  by  about  30%.  The  order  of  the  final  state  space  system  is  now 
about  three  times  as  laxge  as  the  basic  model  due  to  the  two  additional  sets  of  differential 
equations  (typically  60-80  states).  As  before,  only  four  modes  are  needed  to  capture  the 
dynamics  up  to  three  rotor  revolutions  so  that  the  model  can  again  be  reduced  to  eight 
states,  but  this  time  significantly  more  computational  work  is  needed  before  we  can  get 
down  to  this  number.  Solving  for  the  steady  flow  of  the  extended  model  and  computing  the 
state  space  matrices  takes  about  5  times  as  long  as  with  the  basic  model. 
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Table  3.1:  Extended  Hynes-Greitzer  model. 


The  steady  flow  <j>s  and  steady  losses  and  las  are  given  by  the  solution  of 


Pt'1  -  Ps6)  =  +Its  +  ks  + 

(3.65) 

%  /  .  i 

Tra0_  r  +  ZrS' 

(3.66) 

The  complete  state  space  description  is  given  by 

5x  =  A5x  +  B5 7  +  E5 7  4-  G[<5p[l\  ^Pg65]2', 

(3.67) 

_  qp 

<5x  =  [<ty(3),  «r,  <5ls,  <5p£5)]  , 

(3.68) 

A*>  =  Mc-f  -  -d-21^0  ~  Pj,-  ^54 

(3.69) 

Cn 

II 

7^ 

0 

0 

(3.70) 

(3.71) 

-Ay  -A-y  -ty** 

F(4>,%£)I(tM  -\F{4,,)/IjA)  +  D,\  0  0 

F{<f>^)/{r{bs)  0  -[F{<f>s)/(Tibs)  +  Dg}  0 

A^fj,  0  0  —  l/(^a^t^t^so) 

(3.72) 
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Table  3.2:  Three-stage  compressor  model  parameters. 


Nominal  operating  speed  uw  =  72  m/s  at  2400  revolutions  per  minute. 
Geometric  parameters. 


=  0.6792 

fj.s  =  0.3335 

/xa  =  0.2863 

Hi  =  0.0709 

bt  =  0.1184 

bs  =  0.1079 

l\  =  2.9923 

k  =  1-5294 

bG  =  0.1631 

CO 

o 

1 

II 

£ 

H 

f  =  286  mm 

Pressure  rise  characteristics  (measured  by  Haynes  [20]);  see  Figure  3.2  p.67. 
0  =  -10.07 02  +  9.43060  -  1.1849 
0i  =  — 15.534103  +  24.123802  -  15.02620  +  4.6951 

^  =  2.888O02  -  3.65500  +  0.8251 
d'y 

Loss  parameters. 

r  =  0.75 
Tf  =  1.5 


(3.74) 


3.5  Complex  Transfer  Functions 

The  state  space  model  derived  in  the  previous  section  tells  us  how  the  system  behaves  in 
the  time  domain.  It  is  often  convenient  to  analyze  the  system  in  the  frequency  domain 
by  looking  at  transfer  functions.  For  rotating  systems  it  is  attractive  to  view  the  inputs 
and  outputs  as  phasors,  and  phasors  are  best  represented  as  complex  numbers.  In  this 
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Figure  3.2:  Three-stage  compressor  pressure  rise  characteristics. 

section  we  show  how  to  transform  the  real-input  real-output  state  space  description  to  a 
complex-input  complex-output  description.  The  resulting  transfer  function  will  be  called  a 
complex  transfer  function.  This  is  a  slight  misnomer  because  a  transfer  function  is  always 
complex  —  the  term  complex  refers  to  the  fact  that  we  are  using  complex  input  and  output 
signals. 

The  nth  complex  Fourier  coefficient  of  the  flow  perturbation  is  defined  by 

S4>n{r)  =  8<j>cn  ij )  +  i6<j>sn(r).  (3.75) 

Coefficients  for  negative  values  of  n  are  conjugates  of  the  corresponding  positive  valued 
coefficients.  We  note  that  this  definition  differs  from  the  usual  relation  between  real  and 
complex  Fourier  coefficients.  This  nonstandard  definition  is  used  so  that  the  resulting  phasor 
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rotates  in  the  same  direction  as  the  rotor.  In  exactly  the  same  way  we  define 

<57n  =  <57cn  +  i57sn  (3.76) 

for  the  AGV  deflections.  The  expressions  for  computing  the  transfer  function  from  Sjn 
to  6(f>n  given  the  real  input-output  state  description  is  derived  in  Appendix  C.l.  All  the 
transfer  functions  shown  will  be  for  complex  inputs  and  outputs. 

Power  density  spectra  (PSDs)  are  used  quite  often  to  analyze  the  experimental  data  and 
here,  too,  we  will  compute  the  temporal  spectral  densities  of  the  complex  coefficients  defined 
in  Equation  (3.75).  The  Fourier  transform  of  a  complex  function  does  not  have  conjugate 
symmetry  so  that  the  PSD  at  positive  frequencies  will  be  different  from  the  corresponding 
negative  frequencies.  See  also  Paduano  [37]. 


3.6  Modal  Analysis 

In  this  section  we  discuss  the  distorted  flow  mode  shapes  and  compare  them  to  the  uniform 
flow  mode  shapes. 

The  pole  pattern  of  the  neutrally  stable  system  as  predicted  by  the  extended  model 
is  shown  in  Figure  3.3.  The  corresponding  mode  number  is  shown  next  to  the  eigenvalues. 
We  now  explain  how  the  modes  were  numbered  by  considering  the  mode  shapes. 

The  mode  shapes  and  relative  strength  of  the  harmonics  in  the  four  dominant  modes  are 
shown  in  Figure  3.4.  The  dominant  modes  axe  defined  as  those  with  eigenvalues  closest  to  the 
origin.  The  mode  shapes  shown  in  these  figures  were  all  computed  from  the  eigenvectors  of 
the  A— matrix  (see  Chapter  3).  This  choice  was  made  as  the  eigenvectors  are  time  invariant 
while  the  mode  shapes  axe  time  varying.  Putting  it  another  way,  the  mode  shapes  obtained 
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Real  A 


Figure  3.3:  Eigenvalues  of  marginally  stable  system,  1.9  dynamic  head  distortion. 

The  numbers  0,  1,  2,  and  3  to  the  left  of  the  eigenvalues  indicate  the  mode 
number.  Note  that  only  the  eigenvalues  close  »-o  the  origin  are  shown  in  this 
figure. 

from  the  eigenvectors  correspond  to  the  mode  shapes  at  time  zero.  In  these  figures  the 
magnitude  of  the  largest  harmonic  (or  plenum  pressure  perturbation  <5p|5))  was  normalized. 

The  mode  with  eigenvalue  A  =  -0.31  +  0.46i  (top  right  hand  graph)  has  a  strong  first 
harmonic  component.  However,  the  large  contribution  of  the  plenum  pressure  perturbation 
Spi5)  and  presence  of  the  zeroth  harmonic  suggest  that  this  is  a  surge-like  mode,  and  we  will 
therefore  refer  to  it  as  the  zeroth  mode.  Note  that  the  second  harmonic  also  has  significant 
presence  in  the  mode.  The  real  and  imaginary  parts  of  the  corresponding  mode  shape  at 
time  r  =  0  is  shown  in  the  top  left  graph. 
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Figure  3.4:  Mode  shapes  (left)  and  corresponding  harmonic  contribution  (right). 
The  solid  and  dashed  lines  show  the  real  and  imaginary  parts  of  the  mode 


shape  respectively.  The  harmonic  number  indicated  by  p  is  the  plenum  pressure 
perturbations  5pi*\  1.9  dynamic  head  distortion. 


The  neutrally  stable  mode  shown  in  the  second  set  of  graphs  from  the  top  in  Figure  3.4 
also  has  strong  first  harmonic  content  but  the  magnitude  of  the  plenum  pressure  perturba¬ 
tion  5pi5)  is  smaller.  We  thus  call  this  the  first  mode.  The  zeroth  and  second  harmonics 
also  have  substantial  contributions  to  this  mode. 

For  the  other  two  modes  shown  at  the  bottom  of  the  figure  the  contributions  of  the 
zeroth  harmonic  and  5pi$)  decreased  greatly  and  the  modes  consist  mainly  of  the  dominant 
harmonic  plus  the  adjacent  harmonics.  These  two  modes  will  be  called  the  second  and 
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third  modes  respectively.  The  higher  spatial  frequency  content  is  clearly  visible  in  the 
mode  shapes  of  the  second  and  third  harmonics. 

In  all  four  cases  the  mode  shapes  are  not  purely  sinusoidal  and  have  large  deviations 
outside  the  distorted  region  120°  <9  <  240°.  In  Chapter  4  we  saw  that  the  wave  envelope 
was  also  large  outside  the  distorted  region.  With  distortion,  irrespective  of  the  mode  that 
is  excited,  the  envelope  of  the  small  perturbations  will  be  nonuniform  around  the  annulus. 
We  note  that  the  real  and  imaginary  parts  reflect  the  mode  shape  at  two  different  points  in 
time  and  thus  give  an  idea  of  how  the  mode  shapes  change  with  time. 

Neither  the  shapes  nor  the  harmonic  content  of  the  modes  change  significantly  as  the 
extent  of  the  distortion  is  changed  to  180°.  The  magnitude  of  the  distortion  has  a  much 
larger  effect.  Figures  3.5  and  3.6  show  the  modes  and  harmonic  contribution  for  0.8  dynamic 
head  distortion  and  for  uniform  flow  respectively.  The  order  of  the  modes  in  these  figures 
is  the  same  as  in  Figure  3.4.  For  the  0.8  dynamic  head  distortion  the  contribution  of  the 
first  harmonic  to  the  zeroth  mode  is  less  than  half  that  with  the  1.9  dynamic  head  distortion 
because  the  smaller  magnitude  distortion  introduces  less  coupling  between  the  harmonics. 
Similar  arguments  hold  for  the  other  modes.  The  mode  shapes  are  also  more  sinusoidal 
compared  to  1.9  dynamic  head  distortion.  For  uniform  flow  the  mode  shapes  are  purely 
sinusoidal  with  no  coupling  between  the  harmonics. 


3.7  Why  is  Distortion  a  Hard  Problem? 

In  this  section  we  describe  the  main  reasons  that  make  control  of  rotating  stall  a  harder 
problem  in  the  presence  of  distortion. 

Figure  3.7  shows  the  variation  of  the  real  part  of  the  most  unstable  pole  as  a  function 
of  the  distortion  extent.  The  flow  coefficient  was  held  constant  at  1%  above  the  uniform 
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Figure  3.5:  Mode  shapes  (left)  and  corresponding  harmonic  contribution  (right). 
The  solid  and  dashed  lines  show  the  real  and  imaginary  parts  of  the  mode 
shape  respectively.  The  harmonic  number  indicated  by  p  is  the  plenum  pressure 
perturbations  Spi*\  0.8  dynamic  head  distortion. 


flow  open  loop  stall  point  while  the  extent  was  varied  from  0°to  360°. 

Without  distortion  (extent=0)  the  system  is  stable,  that  is,  the  real  part  of  the  eigen¬ 
value  A  is  negative.  As  the  extent  of  the  distortion  is  increased  the  pole  moves  closer  to 
the  itjj— axis  and  goes  unstable  at  approximately  70°  and  remains  unstable  up  to  an  extent 
of  230°.  The  real  part  of  the  pole  assumes  its  maximum  positive  (most  unstable)  value 
at  approximately  135°.  This  is  a  worst  case  situation  and  was  the  motivation  for  using  a 
distortion  with  120°  extent. 
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Figure  3.6:  Mode  shapes  (left)  and  corresponding  harmonic  contribution  (right). 
The  solid  and  dashed  lines  show  the  real  and  imaginary  parts  of  the  mode 
shape  respectively.  The  harmonic  number  indicated  by  p  is  the  plenum  pressure 
perturbations  5pi*\  uniform  flow. 


Also  shown  in  the  figure  is  the  real  part  of  the  dominant  zero.  The  zero  shows  a  similar 
trend  but  crosses  the  iu— axis  at  a  smaller  extent  and  stays  in  the  right  half  plane  over 
a  larger  range  of  extents.  Zeros  in  the  right  half  plane  (that  is,  the  real  part  of  the  zero 
is  positive)  are  called  nonminimum  phase  zeros  and  have  an  adverse  effect  on  closed  loop 
system  performance. 

Nonminimum  phase  zeros  limit  the  performance  of  a  feedback  system  in  the  following 
sense.  Let  the  transfer  function  of  the  system  and  controller  be  denoted  by  G{iui)  and 
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Figure  3.7:  Effect  of  distortion  extent  on  the  most  unstable  pole  and  dominant 
zero,  <f)s0  =  0.465,  1.9  dynamic  head  distortion. 

K(iu)  respectively,  and  let  \G{iu)K(iu)\  >  k  >  1  for  0  <  u  <  u>c.  The  gain-bandwidth 
product  kuc  is  then  limited  and  the  frequency  of  the  nonminimum  phase  zero  is  an  upper 
bound  for  uc.  See  Maciejowski  [30].  A  limit  in  the  gain-bandwidth  product  often  translates 
into  a  limit  on  the  bandwidth  of  the  closed  loop  system.  Therefore,  an  increase  in  k  will 
result  in  more  attenuation  of  disturbances,  but  over  a  smaller  range  of  frequencies. 

Systems  with  nonminimum  phase  zeros  often  have  poor  rejection  of  measurement  noise 
[30].  We  showed  in  Section  4.2.2  that  the  magnitude  of  the  noise  was  larger  for  distorted 
flow  and  we  can  expect  that  this  will  have  a  negative  effect  on  controller  performance. 

We  saw  in  Figure  3.7  that  the  dominant  zero  is  nonminimum  phase  for  a  large  range 
of  extents.  Further,  Figure  3.8  shows  the  dominant  zeros  for  uniform  flow,  0.8  and  1.9 
dynamic  head  distortions  (extent=120°).  The  zero  becomes  nonminimum  phase  at  higher 
mass  flows  for  larger  distortions.  The  zeros  of  the  uniform  flow  and  0.8  dynamic  head 
distortion  systems  become  nonminimum  phase  at  approximately  <f>  —  0.464  and  (f>  =  0.467 
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Figure  3.8:  Real  part  of  dominant  zeros  as  function  of  <j>  for  uniform  flow,  0.8 
and  1.9  dynamic  head  distortions. 


respectively.  Over  the  range  of  flow  coefficients  shown^  the  zero  of  the  system  with  a  1.9 
dynamic  head  distortion  is  always  nonminimum  phase. 

Whether  a  zero  is  minimum  or  nonminimum  phase  is  strongly  affected  by  the  type 
and  configuration  of  sensors  and  actuators.  For  example,  the  model  predicts  that  measur¬ 
ing  velocity  between  the  first  rotor  and  stator  instead  of  upstream  of  the  compressor  will 
increase  the  frequency  of  the  nonminimum  phase  zero  by  50%  for  the  1.9  dynamic  head 
distortion.  An  increased  frequency  implies  an  increased  gain-bandwidth  product,  thus  we 
expect  an  increase  in  the  performance  of  the  controller.  Whether  this  results  in  measurable 
improvement  in  performance  is  investigated  in  Chapter  5. 

In  summary,  the  detrimental  effect  of  distortion  is  twofold.  First,  a  compressor  oper¬ 
ating  in  the  presence  of  distortion  becomes  unstable  at  higher  flow  coefficients.  Second,  the 
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dominant  zeros  become  nomninimum  phase  at  higher  flow  coefficients,  limiting  the  closed 
loop  bandwidth  and  increasing  the  sensitivity  to  measurement  noise. 
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4  Experimental  Results:  Open  Loop 


In  this  chapter  we  assess  the  predictive  capability  of  the  extended  Hynes-Greitzer  model 
for  steady  state  and  dynamic  compressor  behavior.  Comparison  is  made  with  compressor 
behavior  with  distortions  of  two  magnitudes,  0.8  and  1.9  dynamic  heads.  Because  the  results 
were  generally  similar,  unless  stated  otherwise,  all  the  discussions  in  this  chapter  is  for  the 
1.9  dynamic  head  distortion. 

The  chapter  is  divided  into  four  parts.  The  first  part  presents  the  steady  state  velocity 
and  pressure  profiles,  and  pressure  rise  characteristics.  Next,  the  unforced  dynamics  of 
the  compression  system  is  analyzed  in  the  time  and  frequency  domain  (spectral  analysis). 
This  is  followed  by  a  comparison  of  the  predicted  and  measured  forced  response  of  the 
compression  system.  The  chapter  ends  with  a  summary. 


4.1  Steady  State  Results 

In  this  section  we  assess  the  steady  state  predictive  capability  of  the  model.  Steady  state 
results  are  important  for  two  reasons.  First,  the  nonlinear  model  is  linearized  about  the 
nonuniform  steady  velocity  so  that  an  accurate  estimate  of  the  velocity  is  desirable.  Second, 
the  steady  pressure  rise  in  the  presence  of  distortion  is  an  important  measure  of  performance. 


4.1.1  Steady  state  profiles 

The  steady  flow  obtained  by  solving  the  nonlinear  equations  (3.65)  and  (3.66)  gives  the  flow 
at  the  inlet  of  the  first  rotor.  Pressure  measurements,  however,  were  carried  out  at  a  station 
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approximately  0.6  mean  radii  upstream  of  the  compressor  face  (see  Figure  2.4).  We  therefore 
comment  first  on  the  relation  between  these  measurements  and  the  flow  at  the  inlet  of  the 
first  rotor.  The  upstream  static  pressure  decays  exponentially.  We  can  compute  the  static 
pressure  and  velocity  at  the  sensor  location  by  multiplying  the  nth  harmonic  of  the  static 
pressure  at  the  inlet  of  the  first  rotor  by  enXm.  Predicted  and  measured  static  pressures  axe 
shown  in  Figure  4.1.  The  experiment  and  analytical  profile  (with  mean  values  subtracted) 
of  the  XHG  model  are  shown.  The  profiles  have  the  same  shape  and  the  magnitudes  agree 
well. 


Gong  [15]  has  relaxed  the  assumption  that  the  upstream  flow  field  can  be  modelled  in 
a  linearized  way  and  used  a  full  nonlinear  Euler  model  instead;  we  will  refer  to  this  model 
as  the  XHG-Euler  model.  The  result  from  his  calculations  is  also  shown  in  Figure  4.1  and 
again  there  is  good  agreement  with  the  experimental  results.  This  shows  that  the  nonlinear 
effects  in  the  unsteady  flow  field  are  small.  In  all  the  figures  the  extended  Hynes-Greitzer 
and  Euler  models  will  be  indicated  by  XHG  and  XHG-Euler  respectively. 

The  predicted  and  measured  steady  state  velocity  profiles  at  the  measurement  location 
x  —  xm  axe  shown  Figure  4.2.  The  velocity  predicted  by  the  XHG-Euler  model  matches 
the  measured  profile  in  magnitude  well.  The  width  of  the  distorted  region  is  slightly  wider 
than  the  predicted  width.  The  model  predicts  a  lower  velocity  in  the  distorted  region. 
This  can  be  explained  by  comparing  the  predicted  and  measured  total  and  static  pressure 
shown  in  Figure  4.3.  In  this  figure  we  see  that  the  mean  value  of  the  measured  static 
pressure  is  lower  than  the  predicted  pressure  by  10%  so  that  the  predicted  difference  between 
total  and  static  pressures  in  the  distorted  region  is  smaller  than  measured,  resulting  in  a 
lower  estimated  velocity.  This  is  especially  true  at  approximately  200°  where  the  predicted 
velocity  (Figure  4.2  shows  a  laxge  drop.  The  lower  static  pressure  may  be  due  to  radial 
redistribution  upstream  of  the  compressor. 
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Figure  4.1:  Static  pressure  at  x  =  xm,</>s0  =  0.5,  1.9  dynamic  head  distortion. 
The  mean  values  have  been  subtracted  to  show  the  variation  around  the  annulus. 
The  distortion  screen  blocked  the  flow  in  the  range  120°  <6.  <  240°  shown  by 
the  dotted  lines. 


In  Figure  4.3  we  also  show  the  deviation  of  the  exit  static  pressure  ps4).  The  assump¬ 
tion  was  made  that  the  exit  static  pressure  is  uniform  around  the  annulus  and  we  see  the 
measurements  support  this  assumption. 

In  Figure  4.4  we  compare  the  velocities  at  the  compressor  face  predicted  by  the  XHG 
and  XHG-Euler  models;  no  velocity  measurements  are  available  at  this  location.  The  mag¬ 
nitudes  of  the  velocities  agree  very  well  and  the  only  difference  between  the  profiles  is  due 
to  the  narrowing  of  the  flow  field. 
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o.i  h  —  —  XHG 


o  Experiment; 


Figure  4.2:  Steady  state  velocity  profiles  at  x  =  xm,  <j)so  =  0.5  1.9  dynamic  head 
distortion. 


o  p 

+  p 


Figure  4.3:  Total  and  static  pressures,  </>s0  =  0.5,  1.9  dynamic  head  distortion. 
The  mean  value  of  the  static  pressure  at  the  compressor  exit  ps4)  has  been 
subtracted  to  reveal  the  variation  around  the  annulus.  The  solid  and  dashed 
lines  are  for  the  extended  Hynes-Greitzer  model. 


-4  -3.5  -3  -2.5  -2  -1.5  -1  -0.5  0  0.5 


Axial  position  [f] 

Figure  4.5:  Axial  velocity  contours  predicted  by  XHG-Euler  model.  The  vertical 
dashed  lines  at  x  =  -0.5  and  x  =  0  locate  the  inlet  guide  vanes  and  compressor 
face  respectively. 
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Velocity  contours  in  the  x  —  6  plane  calculated  by  the  XHG-Euler  model  are  shown  in 
Figure  4.5.  The  narrowing  of  the  low  stagnation  pressure  region  is  clearly  visible  between 
roughly  100°  and  220°.  We  also  note  the  slight  asymmetry  in  the  region  of  the  distortion 
screen  at  i  s  —3.5  suggesting  some  interaction  between  the  compressor  and  distortion 
generator.  (The  XHG-Euler  model  assumes  an  infinitely  long  duct  upstream  of  the  screen 
while  the  experimental  setup  has  a  bellmouth.) 

We  can  also  estimate  tangential  velocity  using  the  following  procedure.  We  assume 
that  the  irrotational  part  of  the  velocity  equals  the  difference  between  the  velocities  far 
upstream  and  at  the  sensor  location;  this  is  shown  in  the  graph  at  the  top  of  Figure  4.6.  This 
gives  enough  information  to  solve  for  the  tangential  velocity  which  is  compared  against  the 
tangential  velocity  obtained  from  the  Euler  model  in  the  graph  at  the  bottom  of  Figure  4.6. 
Outside  the  regions  where  the  axial  velocity  4>s  has  large  gradients  and  the  assumption  thus 
not  likely  to  be  correct,  there  is  good  agreement  between  the  experimentally  derived  and 
predicted  velocities.  This  again,  is  a  comparison  that  shows  that  nonlinear  effects  will  be 
very  localized  as  far  as  the  upstream  flow  field. 


4.1.2  Speed  lines 

Measured  and  predicted  compressor  performance  curves  for  the  0.8  and  1.9  dynamic  head 
total  pressure  distortions  are  shown  in  Figures  4.7  and  4.8  respectively.  For  the  0.8  dynamic 
head  distortion  the  loss  in  pressure  rise  and  the  increase  in  stalling  flow  coefficient  are  small. 
The  loss  in  pressure  rise  is  1.6%  and  the  stalling  flow  coefficient  increased  by  0.4%.  Both 
these  numbers  are  predicted  correctly  by  the  XHG  model.  The  large  distortion  clearly  shows 
the  detrimental  effect  of  distortion  —  the  pressure  rise  dropped  by  7.6%  and  the  stalling 
flow  coefficient  increased  by  4.3%.  The  model  predicted  an  increase  of  2.0%  in  stalling  flow 
coefficient  and  a  drop  in  peak  pressure  rise  of  4.4%.  The  characteristic  that  has  been  used 
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*  Experiment,  uniform  flow 

-  Curve  fit  for  uniform  flow 

o  Experiment,  small  distortion 

—  —  extended  Hynes-Greitzer  model 


_ I _ I _ I _ I _ I _ L _ 

0.42  0.44  0.46  0.48  0.5  0.52 


0.54  0.56  0.58 


Figure  4.7;  Compressor  performance  line  for  0.8  dynamic  head  distortion. 
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Figure  4.8:  Compressor  performance  line  for  1.9  dynamic  head  distortion 


was  measured  by  Haynes  [20]  (see  3.2  on  page  66)  over  the  range  0.42  <  <p  <  0.7  and  it 
has  been  assumed  that  the  parabola  that  fits  the  data  over  this  range  of  flow  coefficients  is 
also  valid  at  flow  coefficients  down  to  0.37,  see  Figure  4.4.  The  lower  predicted  stalling  flow 
coefficient  and  smaller  predicted  drop  in  peak  pressure  rise  suggest  that  the  slope  of  the 
characteristic  to  the  left  of  the  peak  (0SO  =  0.460)  may  be  steeper  than  that  of  the  parabola. 
Currently  no  method  exists  to  measure  the  characteristic  over  the  complete  range  of  flow 
coefficients  needed  for  the  large  distortion. 

Table  4.1:  Open  loop  stalling  flow  coefficients.  Percentages  are  computed  rela¬ 
tive  to  (j>  =  0.460.  _ 


Uniform  flow 

120°,  0.8 

120° ,  1.9 

^ataii  XHG  model 
Natali  experiment 

0.460 

0.460 

0.462  0.4% 

0.462  0.4% 

0.469  2.0% 

0.480  4.3% 

The  extended  model  was  used  to  obtain  the  predictions  shown.  The  basic  model 
predicts  the  trends  but  gives  a  higher  stalling  flow  coefficient  than  the  experiment.  Despite 
this  limitation  the  basic  model  will  be  shown  to  be  a  useful  tool  for  designing  controllers  to 
stabilize  the  flow. 


4.2  Dynamics 

4.2.1  Time  Domain  Analysis 

In  this  section  we  compare  the  dynamic  behavior  of  the  compression  system  with  that 
predicted  by  the  XHG  model  using  time  domain  data.  For  the  data  available  the  signal  to 
noise  ratios  (SNRs)  of  the  velocity  measurements  are  less  than  one  so  it  is  often  difficult  to 
distinguish  the  signal  from  the  noise.  The  goal  therefore  is  not  so  much  to  find  how  well  the 
model  predicts  the  behavior  of  the  compression  system  but  to  observe  trends  of  the  small 
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Real  A 


Figure  4.9:  Eigenvalues  of  marginally  stable  system,  1.9  dynamic  head  distortion. 
The  numbers  0,  1,  2,  and  3  to  the  left  of  the  eigenvalues  indicate  the  mode 
number.  Note  that  only  the  eigenvalues  clcce  to  the  origin  are  shown  in  this 
figure. 


perturbation  dynamics  in  the  presence  of  distortion. 

The  pole  pattern  of  the  neutrally  stable  system  as  predicted  by  the  extended  Hynes- 
Greitzer  model  is  shown  in  Figure  4.9.  For  reference  we  list  the  eigenvalues  and  correspond¬ 
ing  mode  number  shown  in  the  figure. 

A0  =  -0.31  +  0.46i  Ai  =  -0.00  +  0.38i 

A2  =  — 0.20  +  0.85z  A3  =  — 0.39  +  1.35t. 

In  this  figure  we  show  only  eigenvalues  with  magnitude  less  than  1.5  uT  but  a  complete  set 
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of  eigenvalues  will  be  shown  later.  We  will  refer  to  the  neutrally  stable  mode  with  natural 
frequency  =  0.38  as  the  first  mode,  similarly  for  others.  The  numbering  of  the  modes 
was  explained  in  Section  3.6. 

We  will  compare  the  behavior  of  the  experimental  system  just  prior  to  stall  with  that 
of  the  unforced  model  when  only  the  first  mode  has  nonzero  initial  conditions.  All  the 
modes  of  the  experimental  system  will  be  excited  by  noise.  The  marginally  stable  mode 
should  dominate  the  behavior  of  the  system,  so  it  should  correspond  to  the  simulation.  The 
homogeneous  system 

Sx  =  ASx  (4.1) 

has  the  solution 

Sx  =  eArSx0.  (4.2) 

We  will  choose  the  initial  state  Sx0  so  that  only  the  first  mode  is  excited.  The  solution  can 
be  simplified  if  we  express  the  A-matrix  in  terms  of  the  eigenvalues  and  eigenvectors.  Let 

A  =  VAV~\  (4.3) 

where  V  =  [vi,  u2, . . .]  is  the  matrix  of  eigenvectors,  (4.4) 

and  A  =  diag[Ai,  A2, ...]  is  the  matrix  of  eigenvalues.  (4.5) 

We  assume  without  loss  of  generality  that  the  first  two  eigenvalues  and  eigenvectors  corre¬ 
spond  to  the  neutrally  stable  mode  so  that  Ai  =  iu>i ,  and  A2  =  With  this  choice 

Sx  =  SR(vi)  cos(wit)  -  Q'(ui)  sin(wir)  (4.6) 

where  SR  and  Q  indicate  the  real  and  imaginary  part  respectively.  Let  the  eigenvector  be 


v\  =  [wo,  vci,vsi,vC2,  vs2, . . .  ]T  then  the  reconstructed  perturbation  is 
6<j>(6 ,  r)  =  3t(t>o)  cos(a;ir)  —  S( vq )  sin^r) 

+  ^2  [3®(vCT»)  cos(wit)  -  S(uCn)  sin(o/’iT )]  cos(n0) 

n>0 

+  ^2  t^(usn)  cos(wir)  -  S(vsn)  sin(wir)]  sin(n0)  (4.7) 

n>  0 

=  5<j>o  +  ^2  ^cn  cos(n^)  +  <^sn  sin(n0).  (4.8) 

n>0 

Notice  that  the  frequency  of  oscillation  of  the  coefficients  5<f> o,  8<f>cn,  and  5<fisn  of  the  nth 
harmonic  is  the  same  for  all  the  harmonics,  see  Equations  (4.7)  and  (4.8).  The  expressions 
for  the  coefficients  can  be  simplified  by  using  standard  trigonometric  identities  to  obtain 

t )  =  aco  cos(wxr)  +  aso  sin(wir) 

+  ^2  acn  cos(n^  —  wiT)  +  asn  sin(n0  -  0>ir) 
n>  0 

+  ^>2  b cn  cos(n&  +  u>iT )  +  bsn  sin(n0  +  u}\t)  (4.9) 

n>  o 

and  the  amplitudes  ac0,  as0,  acn,  asn,  6cn,  and  bsn  depend  the  elements  of  the  eigenvector; 
the  exact  expressions  are  not  relevant  to  our  discussion.  The  important  point  is  that  a  mode 
can  be  decomposed  into  two  waves:  one  travelling  in  the  direction  of  rotor  rotation  and  one 
travelling  in  the  opposite  direction  and,  in  general,  these  waves  have  different  amplitudes 
and  phases.  Because  the  amplitudes  of  the  two  waves  differ  the  wave  shape  will  change  with 
time.  No  assumptions  were  made  about  the  distortion  so  Equations  (4.8)  and  (4.9)  hold  in 
general. 

An  example  of  the  waves  predicted  by  Equation  (4.8)  (or  equivalently  Equation  (4.9)) 
is  shown  in  Figure  4.10.  Five  wave  forms  approximately  equally  spaced  in  time  over  one 
period  are  shown  together  with  the  envelope  of  the  wave.  The  variation  in  magnitude  of 
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the  envelope  around  the  annulus  can  be  explained  by  looking  at  the  slope  shown  in  Fig¬ 
ure  4.11.  The  magnitude  of  the  envelope  decreases  when  the  slope  is  negative  and  increases 
when  the  slope  is  positive.  Therefore,  the  magnitude  of  the  envelope  has  a  maximum  where 
the  slope  changes  from  positive  to  negative. 

A  typical  data  set  of  stall  inception  is  shown  in  Figure  4.12. 

The  origin  of  the  time  axis  has  been  chosen  as  the  point  where  the  velocities  exceed 
an  arbitrary  large  limit.  The  exact  value  of  the  origin  is  not  important  as  we  are  interested 
in  the  period  prior  to  stall.  Two  properties  predicted  by  the  model  are  presented  in  the 
figure.  Two  parallel  lines  drawn  through  the  peaks  of  the  perturbations  are  used  to  help 
show  that  the  wave  is  travelling  around  the  annulus.  An  estimate  of  the  speed  at  which  the 
wave  rotates  can  be  obtained  by  counting  the  number  of  peaks  in  the  trace  at  the  bottom 
of  the  figure;  between  r  =  —20  and  r  =  — 10  there  are  a  little  more  than  four  peaks  giving 
a  speed  of  just  over  0.4  wr  which  agrees  well  with  the  predicted  0.38  lut. 

Second,  if  we  look  around  the  annulus  at  a  fixed  time  (prior  to  stall)  we  see  that  the 
perturbation  magnitude  is  large  for  small  values  of  6,  decreases  as  9  increases,  is  small  for 
120°  <  6  <  240°,  and  increases  again  towards  360°;  this  is  predicted  by  the  model,  see 
Figure  4.10. 

Analysis  of  the  perturbations  around  the  annulus  has  shown  that  they  are  not  dis¬ 
tributed  symmetrically  about  the  mean  values.  The  standard  deviations  (STDs)  for  positive 
and  negative  values  of  the  perturbations  are  shown  in  Figure  4.13  where  we  see  that  for 
negative  values  of  the  perturbations  the  STDs  match  the  envelope  well,  but  for  positive 
values  of  the  perturbations  the  STD  shows  little  variation  around  the  annulus.  This  asym¬ 
metric  behavior  is  not  predicted  by  the  model  and  shows  either  that  there  are  dynamics 
not  captured  by  the  model  or  nonlinear  behavior.  For  insight  into  the  prestall  behavior 
we  decompose  the  perturbations  64>{9,t)  into  a  Fourier  series  with  complex  harmonics  as 
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Figure  4.12:  Circumferential 


Figure  4.13:  Standard  deviation  of  the  perturbations  around  annulus  for  S(f>  >  0 
and  5(j>  <  0.  The  standard  deviations  obtained  from  the  different  hot-wires  were 
connected  by  straight  lines  to  show  the  trend. 


discussed  in  Section  3.5..  The  nth  harmonic  is  indicated  by  54>n{r).  The  magnitudes  of  the 
zeroth,  first,  second,  and  third  harmonics  are  shown  in  the  graph  at  the  top  of  Figure  4.14. 
The  nth  harmonic  has  been  scaled  by  e+nlI,nl  so  the  magnitudes  shown  are  those  seen  by 
the  compressor.  In  this  figure  we  see  that  the  magnitudes  of  the  higher  harmonics  increase 
gradually  from  \5(f>n  |  «  0.035  for  approximately  20  rotor  revolutions  prior  to  stall. 

The  unwrapped  phases  of  the  harmonics  are  also  shown  in  Figure  4.14  where  the 
travelling  of  the  first  harmonic  at  0.43  u>T  is  clearly  visible.  The  phase  of  the  second  harmonic 
is  more  erratic  but  shows  a  definitive  trend;  the  corresponding  velocity  of  0.20  ljt  is  slightly 
lower  than  0.43/2  =  0.215.  The  phase  speeds  prior  to  stall  were  computed  with  least  squares 
fits  over  the  range  —15  <  r  <  —5.  Poorer  estimates  for  the  higher  harmonics  must  be 
expected  because  the  higher  harmonics  decay  exponentially  and  thus  have  lower  SNRs. 
Phase  unwrapping  is  sensitive  to  noise  and  so  the  velocities  obtained  from  the  phases  should 
be  considered  rough  estimates.  Even  so,  the  velocities  for  the  first  and  second  harmonics 
agree  well  with  the  predicted  values  of  0.38  uT  and  0.19  uT  respectively. 


4.2.2  Spectral  Analysis 

The  power  density  spectrum  (PSD)  is  a  powerful  tool  for  analyzing  the  behavior  of  the 
system.  The  PSDs  used  here  were  computed  using  Welch’s  method  [36]  with  a  Hanning 
window. 

Figure  4.15  show  the  PSDs  of  the  complex  harmonics  S<po-5(pz-  Except  for  the  third 
harmonic  the  first  mode  shows  up  strongly  at  u>i  =  0.44 uT  in  the  PSDs  indicating  there  is  a 
wave  travelling  around  the  annulus  at  44%  of  the  rotor  speed.  The  travelling  wave  and  peak 
in  the  PSDs  of  the  harmonics  are  predicted  by  the  model  (see  Equations  (4.8)  and  (4.9)) 
and  shows  the  strong  coupling  between  the  different  harmonics,  including  the  zeroth.  The 
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coupling  results  from  the  nonuniform  steady  flow  and  does  not  occur  with  uniform  mean 
flow;  for  uniform  only  the  nth  harmonic  is  present  in  the  nth  mode.  The  PSDs  for  uniform 
flow  are  shown  in  Figure  4.16.  The  first  mode  at  0.28wr  is  clearly  visible  in  the  PSD  for 
6(j>i  and  is  barely  visible  in  the  PSDs  of  the  other  harmonics. 

A  comparison  of  the  noise  floors  for  distorted  flow  (Figure  4.15)  with  the  corresponding 
uniform  flow  (Figure  4.16)  shows  that  the  magnitude  of  the  noise  is  bigger  for  distorted 
flow.  If  the  peak  at  1  u>T  is  excluded  from  the  calculations,  the  distorted  flow  SNR  for 
the  zeroth,  first,  and  second  harmonics  are  0.34(— 9.4  dB),  0.63(— 4  dB),  and  0.29(— 10  dB) 
respectively.  The  SNR  for  the  first  harmonic  for  uniform  flow  is  1.25(1.9  dB)  that  is  better 
than  the  corresponding  distorted  case.  The  reason  for  the  lower  SNR  when  distortion  is 
present  is  not  known.  The  PSDs  were  computed  over  a  period  of  2400  rotor  revolutions 
prior  to  stall.  Directly  prior  to  stall  the  magnitudes  of  the  perturbations  tend  to  increase 
so  that  the  SNR  improves  and  it  becomes  possible  to  distinguish  the  signal  from  the  noise 
without  additional  signal  processing  as  we  have  seen  in  Figure  4.12. 

In  Figure  4.15  we  also  see  a  large  peak  at  exactly  lwr.  Comparison  with  the  uniform 
flow  PSDs  in  Figure  4.15  shows  that  the  lwr  peak  increased  by  approximately  20,  10,  40, 
and  30  dB  for  the  zeroth  to  third  harmonics  respectively  when  distortion  was  added.  It 
thus  appears  that  this  disturbance  is  related  to  circumferential  nonuniform  flow.  We  note 
that  the  PSDs  also  have  peaks  at  2  and  3  ujt  which  are  not  shown. 

In  Figure  4.15  we  also  see  peaks  at  0.56 ujt  in  the  graphs  for  5(f) o,  5d>\ ,  and  5<f>2  and 
at  1.44  uT  for  5<f> i,  5<f> 2,  and  dfa.  The  frequencies  0.56  ujt  and  1.44  u>r  are  exactly  u>r  ±  u>i 
and  shows  coupling  between  the  first  mode  and  the  disturbance  at  one  rotor  revolution. 
This  coupling  has  the  form  of  modulation  (or  multiplication)  in  the  time  domain,  resulting 
in  a  sum  and  difference  frequency  in  the  frequency  domain.  The  origin  of  the  one  rotor 
revolution  mode  and  modulation  is  not  known. 
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We  also  note  that  the  PSDs  at  negative  frequencies  (shown  in  dashed  lines)  differ  from 
those  at  the  corresponding  positive  frequencies.  Generally  there  is  very  little  activity  at 
negative  frequencies. 


4.3  Process  Noise  Model 

In  the  previous  section  we  saw  that  there  was  substantial  noise  present  in  the  compression 
system.  Two  external  disturbances  forcing  the  system,  (far  upstream)  and  <5ps!)  (at 
throttle  exit),  were  included  in  the  model  (see  Chapter  3)  as  process  noise.  These  noise 
sources,  however,  do  not  describe  the  noise  seen  in  the  measurements. 

Figure  4.17  shows  the  analytical  and  measured  PSD  of  the  first  harmonic  with  5p[1] 
as  the  process  noise  for  the  1.9  dynamic  head  distortion.  A  logarithmic  frequency  axis  was 
chosen  to  show  the  trend.  Except  for  the  first  mode  peak,  the  magnitude  of  the  PSD  has  a 
slope  of  20  dB/decade  as  indicated  by  the  dashed  line.  The  analytical  PSD  has  a  constant 
slope  at  frequencies  below  0.1  ur  and  rolls  off  very  fast  for  frequencies  above  that  of  the  first 
mode.  For  the  analytical  PSD  it  was  assumed  the  noise  6p[l)  is  white. 

The  20  dB/decade  slope  suggests  that  the  integral  of  white  noise  must  be  used  as 
process  noise.  However,  this  will  not  give  a  slope  of  20  dB/decade,  but  40  dB/decade 
because  the  square  of  the  magnitude  of  the  transfer  function  is  used  in  the  computation 
of  the  PSD.  The  only  way  to  obtain  a  20  dB/decade  slope  at  low  frequencies  is  to  assume 
that  the  PSD  of  the  noise  has  a  1/u  slope.  If  this  is  the  case  the  discrepancy  between  the 
measured  and  analytical  PSDs  at  frequencies  above  that  of  the  first  mode  will  be  larger. 

It  may  be  possible  that  the  noise  seen  at  low  frequencies  is  due  to  process  noise  while 
that  at  frequencies  above  uj\  is  measurement  noise  (not  included  in  the  model).  Both  these 
noises  will  have  to  have  1  jm  slope.  Whether  this  explanation  holds  is  not  known. 
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Figure  4.17:  Experimental  and  analytical  PSDs,  1.9  dynamic  head  distortion. 


Similar  results  hold  for  the  other  harmonics  and  process  noise  p(s6)  and,  in  summary,  the 
process  noise  can  at  most  be  in  part  responsible  for  the  noise  seen  in  the  measured  PSDs. 


4.4  Transfer  Functions 

The  state  space  description  and  corresponding  time  domain  analysis  is  one  way  of  looking 
at  signals.  An  alternative  way  of  looking  at  the  system  is  in  the  frequency  domain.  Let  the 
state  space  description  be  given  by 


6x  =  A6x  +  B8i  (4.10) 

8<f>  =  CSx  +  D8~f.  (4.11) 

We  assume  that  the  dependence  on  S'y  has  already  been  accounted  for;  see  discussion 
on  page  59.  Taking  the  Laplace  transforms  of  these  equations  and  substituting  the  first 
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equation  into  the  second  we  obtain 


6<f>(iu))  =  [C{iul  -  A)~lB  +  D]  (4.12) 

=  [CV(iuI  -  A 4-  D\  S^(iu})  (from  Equation  (4.3))  (4.13) 

=  G(iu)  (4-14) 

The  matrix  G(iu )  is  known  as  the  transfer  function  matrix.  The  transfer  function  can  be 
written  in  terms  of  the  eigenvalues  A  and  the  eigenvectors  (or  mode  shapes)  V.  Therefore,  if 
there  is  good  match  between  the  predicted  and  measured  transfer  functions  we  can  conclude 
there  is  good  match  between  the  predicted  and  measured  mode  shapes. 

The  method  used  for  measuring  the  transfer  function  is  discussed  in  Appendix  C.2.  In 
that  appendix  we  also  discuss  the  excellent  noise  rejection  and  robustness  against  nonlinear¬ 
ities  of  the  measuring  technique,  with  the  conclusion  that  examination  of  transfer  functions 
is  the  most  reliable  way  to  determine  the  accuracy  of  the  model. 


All  the  transfer  functions  shown  will  be  for  complex  inputs  and  complex  outputs  as 
discussed  in  Section  3.5  and  Appendix  C.l. 

We  will  show  transfer  functions  from  the  first  three  harmonics  of  £7  to  the  first  three 


harmonics  of  Sep  so  that  G  is  a  square 
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Magnitude  [dB] 


57o  — ►  6<t> 0  ^7i  ->  <5^0 


Frequency  [<jr]  Frequency  [wr] 


Figure  4.18:  Transfer  functions  <f>  =  0.500,  1.9  dynamic  head  distortion.  Solid 
line=predicted,  dots=experiment,  dashed  line=uniform  flow  experiment. 
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Phase  [°] 


J72  -4 


($73  -4  <ty0 


Frequency  [u;r] 


Frequency  [ur\ 


Figure  4.19:  Transfer  functions  <f>  —  0.500,  1.9  dynamic  head  distortion.  Solid 
line=predicted,  dots=experiment,  dashed  line=uniform  flow  experiment. 
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<$70  &<t> 2 


6*yi  —y5(j>  2 


^70  ->  ^3 


#7i  — >-  S(j> 3 


Figure  4.20:  Transfer  functions  <f)  =  0.500,  1.9  dynamic  head  distortion.  Solid 
I  ine= predicted,  dots=experiment,  dashed  line=uniform  flow  experiment. 
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Phase  [°] 


Frequency  [wr]  Frequency  [cur] 

Figure  4.21:  Transfer  functions  <p  =  0.500,  1.9  dynamic  head  distortion.  Solid 
line=predicted,  dots=experiment,  dashed  line=uniform  flow  experiment. 
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Figures  4.18  -  4.21  show  the  set  of  transfer  functions  for  the  1.9  dynamic  head  distortion 
at  <fis0  =  0.500.  This  is  approximately  17%  above  the  stalling  flow  coefficient  so  that  we 
expect  the  system  to  exhibit  behavior  that  is  well  damped.  In  these  figures  the  transfer 
function  gpq  from  6jq  to  8<pp  is  indicated  by  8~fq  — >•  8<pp.  Four  transfer  functions  are  shown 
on  each  page.  The  magnitude  of  each  transfer  function  is  shown  with  the  corresponding 
phase  directly  below  it.  The  transfer  functions  along  the  diagonal  of  G  also  include  the 
measured  uniform  flow  transfer  functions  (shown  with  dashed  lines);  for  uniform  flow  all 
off-diagonal  transfer  functions  are  zero. 

The  first  harmonic  has  a  strong  presence  in  the  first  mode  so  we  will  concentrate 
the  discussion  on  the  transfer  functions  to  and  from  the  first  harmonic.  Peaks  in  the 
transfer  function  magnitudes  indicate  lightly  damped  modes  at  which  the  system  will  tend  to 
resonate.  For  example,  the  measured  transfer  function  from  £71  — >■  8<f>  1  at  the  bottom  right 
hand  corner  of  Figure  4.18  has  a  small  peak  at  approximately  0.38  wr  while  the  predicted 
peak  is  at  0.40  wr.  The  model  thus  predicted  the  frequency  of  the  first  mode  with  an  error 
of  only  0.02  u>r,  that  is,  to  within  2%  of  the  rotor  frequency. 

The  peak  at  u>\  =  0.38  tur  is  also  visible  in  the  magnitudes  of  the  transfer  functions 
goo  and  522  indicating  that  these  harmonics  are  also  present  in  the  first  mode.  The  smaller 
peaks  in  goo  and  322  can  be  explained  by  looking  at  the  predicted  pole-zero  maps  of  the 
individual  transfer  functions  shown  in  Figure  4.22. 

For  goo  we  see  that  there  is  a  zero  close  to  the  first  mode  pole  so  the  transfer  function 
magnitude  will  not  show  a  large  peak  in  the  corresponding  frequency  range.  For  522  there  is 
an  almost  exact  pole-zero  cancellation,  while  <733  has  an  exact  pole-zero  cancellation  with  no 
peak  in  the  transfer  function  magnitude.  We  also  see  in  these  pole-zero  maps  that  the  poles 
at  negative  frequencies  are  cancelled,  or  almost  cancelled,  by  zeros  and  thus  the  transfer 
functions  at  negative  frequencies  do  not  have  any  large  resonant  peaks  and  are  therefore  not 
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Imag  A  [wr]  j  tImag  A  [u>r] 


Real  A  [wr]  Real  A  [wr] 


-1  -0.75  -0.5  -0.25  0  -1  -0.75  -0.5  -0.25  0 

Real  A  [wr]  Real  A  [wr] 


Figure  4.22:  Pole-zero  maps  of  gjj,j  =  0,1, 2, 3.  <j>  —  0.500,  1.9  dynamic  head 
distortion.  x=poles,  o=zeros 
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shown.  The  lack  of  strong  resonance  at  negative  frequencies  was  also  visible  in  the  PSDs 
showed  earlier  on,  see  Figure  4.15. 

Because  of  the  low  SNR,  transfer  functions  with  magnitudes  below  approximately 
—40  dB  will  not  be  very  accurate. 

Comparing  the  uniform  and  distorted  flow  transfer  functions  goo ,  gn,  922,  and  <733  we 
see  that,  except  for  gn,  the  addition  of  distortion  did  not  change  these  transfer  functions 
significantly.  The  magnitude  of  511  for  distorted  flow  is  10  dB  larger  than  the  corresponding 
uniform  flow  transfer  function,  indicating  that  the  eigenvalue  of  the  dominant  mode  has 
moved  closer  to  the  iu— axis,  that  is,  distortion  has  a  destabilizing  effect. 

Valleys  in  the  transfer  functions  indicate  zeros  close  to  the  to;— axis.  Again  this  is 
visible  in  the  transfer  function  gn •  In  this  case  it  appears  from  the  magnitude  of  gn  that 
the  predicted  frequency  of  the  zero  is  lower  than  the  measured  frequency.  The  minimum 
values  of  the  measured  and  predicted  transfer  function  magnitudes  are  at  0.93  uT  and  1.03  u>r 
respectively,  an  error  of  0.1  wr.  The  measured  and  predicted  phases  of  gn  changes  abruptly 
to  within  0.05  wr,  and  we  conclude  that  the  frequency  of  the  zero  is  predicted  accurately  by 
the  model. 

The  predicted  zero  is  minimum  phase,  that  is,  the  real  part  of  the  zero  is  negative  (see 
Figure  4.22)  so  that  the  phase  increases  for  frequencies  close  to  that  of  the  zero  frequency, 
while  the  measured  zero  is  nonminimum  phase  (real  part  of  the  zero  is  positive)  with  a 
corresponding  decrease  in  the  phase.  At  this  flow  coefficient  the  zero  is  very  close  to  the 
iu— axis  as  is  evident  by  the  abrupt  change  in  both  the  measured  and  predicted  phase.  The 
predicted  and  measured  zero  sire  just  to  the  left  and  right  of  the  iu— axis  respectively  so 
that  the  absolute  distance  between  the  zeros  is  small. 

The  magnitudes  of  the  transfer  functions  goi,  gio,  <712,  and  <721  are  in  the  range  —30  dB 
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to  —20  dB  for  frequencies  close  to  that  of  the  first  mode,  while  the  magnitude  of  <711  is 
approximately  —10  dB  at  the  same  frequencies,  indicating  strong  coupling  between  the 
zeroth,  first,  and  second  harmonics.  Except  for  the  transfer  function  from  the  zeroth  to  the 
first  harmonic  there  is  good  agreement  between  the  magnitudes  and  phases  of  the  predicted 
and  measured  transfer  functions.  The  transfer  function  from  the  zeroth  to  the  first  harmonic 
appears  to  have  additional  zeros  not  predicted  by  the  model;  however,  the  general  trend 
of  the  magnitude  and  phase  still  follows  the  measured  transfer  function  well  up  to  lur. 
In  general  it  was  found  that  transfer  functions  from  the  zeroth  harmonic  to  the  higher 
harmonics  were  not  predicted  as  well  as  the  other  transfer  functions;  however,  the  main 
trend  was  always  captured  well. 

We  also  note  in  these  figures  that  transfer  functions  gpq  for  which  \p  —  q\  <  1  are 
predicted  better  than  transfer  functions  for  which  \p  —  q\  >  1.  This  trend  can  be  explained 
by  the  fact  that  the  higher  harmonics  decay  like  enXm  so  that  the  SNR  is  smaller  for  higher 
harmonics  and  it  is  more  difficult  to  measure  these  transfer  functions  accurately.  In  addition, 
for  higher  harmonics  the  approximation  that  the  12  AG  Vs  are  a  continuum  also  becomes 
poorer.  This  too  will  decrease  the  fidelity  of  the  measured  transfer  functions. 

To  determine  if  the  model  captures  the  change  in  dynamic  behavior  of  the  compression 
system  at  different  flow  coefficients,  transfer  functions  were  measured  at  three  different 
operating  points.  A  subset  of  this  data  is  shown  in  Figures  4.23  -  4.26.  The  graphs  at 
the  top,  middle,  and  bottom  show  the  magnitude  (left)  and  phase  (right)  at  0staii  +  17%, 
</>staii  4-  7%,  and  </>stan  —  0.4%  respectively.  In  general  the  trend  was  captured  well  by  the 
model. 

At  the  lowest  flow  coefficient  the  frequency  of  the  first  mode  is  under  predicted  by 
approximately  0.1  ur.  There  are  two  possible  reasons  for  this  discrepancy.  First,  at  <£staii  — 
0.4%  the  steady  flow  spans  the  range  0.36  <(f>%<  0.55.  The  characteristic  has  been  measured 
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Magnitude  [dB]  Magnitude  [dB]  Magnitude  [dB] 


6-yo  6<j>o 


Figure  4.23:  Magnitudes  (left)  and  phase  (right)  of  transfer  functions  from  <570 
to  5<j)o  at  install  +  17%  (top),  (fistaU  +  7%  (middle),  and  &taii  -  0.4%  (bottom), 
1.9  dynamic  head  distortion.  Solid  line=predicted,  dots=experiment. 
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Frequency  [wr] 


Figure  4.24:  Magnitudes  (left)  and  phase  (right)  of  transfer  functions  from  <571 
to  84> 1  at  (ktan  +  17%  (top),  +  7%  (middle),  and  ^aii  -  0.4%  (bottom), 
1.9  dynamic  head  distortion.  Solid  I ine= predicted,  dots=experiment. 
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Figure  4.25:  Magnitudes  (left)  and  phase  (right)  of  transfer  functions  from  £72 
to  54>2  at  0staii  +  17%  (top),  0staii  +  7%  (middle),  and  </>staii  -0.4%  (bottom), 
1.9  dynamic  head  distortion.  Solid  line=predicted,  dots=experiment. 


Magnitude  [dB]  Magnitude  [d B]  Magnitude  [dB] 


<5  72  $<(>  1 


Figure  4.26:  Magnitudes  (left)  and  phase  (right)  of  transfer  functions  from  Sjz 
to  64>z  at  <£staii  +  17%  (top),  </>stai i  +  7%  (middle),  and  <£stall  -  0.4%  (bottom), 
1.9  dynamic  head  distortion.  Solid  line=predicted,  dots=experiment. 
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only  for  <f)s  >  0.42  and  thus  we  are  operating  on  a  part  of  the  characteristic  that  is  not 
known.  For  the  predicted  transfer  functions  shown  in  the  figure  it  has  been  assumed  that 
the  characteristic  for  flow  coefficients  less  than  0.42  is  parabolic,  and  the  actual  slope  may 
be  different.  Second,  when  the  distortion  generator  was  disassembled  after  completion  of 
the  experiments  it  was  found  that  the  screen  was  packed  with  dirt  and  thus  the  distortion 
intensity  was  higher  than  the  value  measured  initially.  The  frequency  of  the  modes  increase 
with  increasing  distortion  intensity,  so  the  dirty  screen  explains,  at  least  in  part,  the  higher 
measured  frequency  of  the  first  mode. 

At  the  lowest  flow  coefficient  shown,  ^staii  —  0.4%,  the  system  was  open  loop  unstable 
so  it  was  necessary  to  use  active  control  to  stabilize  the  system.  Active  control  will  be 
discussed  in  Chapter  5. 

4.5  Summary 

In  this  chapter  we  compared  the  dynamic  behavior  of  small  perturbations  in  both  the  time 
and  frequency  domain.  All  the  major  trends  are  predicted  by  the  model.  Analysis  of 
the  transfer  functions  showed  the  model  accurately  predicted  the  frequencies  of  the  most 
unstable  pole  and  dominant  zero  to  within  0.02  u>T  and  0.1  wr  respectively. 
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5  Active  Control 


5.1  Introduction 

In  this  chapter  we  present  the  first  experimental  results  of  active  control  of  rotating  stall 
in  the  presence  of  large  total  pressure  inlet  distortion.  Several  questions  will  be  addressed: 
Will  modern  controllers,  whose  designs  axe  based  purely  on  the  model,  stabilize  the  system? 
Are  harmonic  feedback  controllers  designed  for  uniform  flow  (used  by  Paduano  [37]  and 
Haynes  [20])  effective  in  the  presence  of  distortion?  Can  the  coupling  between  the  different 
harmonics  due  to  the  distortion  be  used  to  stabilize  the  compression  system  at  lower  flow 
coefficients?  Does  the  model  predict  the  range  of  flow  coefficients  over  which  the  compression 
system  can  be  stabilized? 

The  figure  of  merit  used  to  determine  the  performance  of  a  controller  will  be  the  lowest 
flow  coefficient  at  which  the  compression  system  can  be  stabilized.  This  is  not  the  only  or 
best  measure  of  performance;  other  factors  like  disturbance  rejection,  sensitivity  to  noise, 
and  robustness  to  system  uncertainty  are  of  practical  importance.  Also,  the  effectiveness  of 
the  controller  against  unknown,  time  vaxying  distortions  and  different  operating  conditions 
are  important  for  real  applications.  However,  for  a  given  compressor,  determining  the 
stalling  flow  coefficient  experimentally  is  convenient,  can  be  done  in  a  consistent  way,  and 
has  been  used  successfully  in  the  past  by  Paduano  [37],  Haynes  [20],  and  Gysling  [19]. 

To  answer  the  questions  posed  earlier  four  controllers  were  tested  experimentally.  We 
start  with  a  discussion  of  the  different  controllers.  This  is  followed  by  a  presentation  of  the 
experimental  results,  an  investigation  into  the  loss  of  stability,  and  a  chapter  summary. 
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5.2  Control  Laws 


5.2.1  LQG  Controller 

The  first  controller  is  a  linear  quadratic  Gaussian  (LQG)  controller.  LQG  controllers  axe 
usually  used  in  optimal  control  applications  where  a  quadratic  cost  function  consisting  of 
weighted  sums  of  the  state  perturbations  and  control  energy  is  minimized.  For  our  applica¬ 
tion  the  LQG  design  procedure  was  not  used  for  its  optimal  properties  but  simply  because  it 
provides  a  simple  way  to  design  controllers  for  multi-input  multi-output  systems.  The  basic 
model  was  used  for  the  design  of  the  LQG  controllers.  Because  the  basic  model  does  not 
include  any  model  of  unsteady  viscous  response  the  predicted  compressor  slope  is  steeper 
than  the  effective  slope;  therefore,  the  slope  of  the  characteristic  at  flow  coefficients  to  the 
left  of  the  peak  of  the  characteristic  was  decreased  until  the  predicted  stalling  flow  coeffi¬ 
cient  matched  the  experimental  valued.  With  this  change  it  was  found  that  the  frequency 
of  the  first  mode  was  also  predicted  correctly. 

The  LQG  controllers  are  designed  in  two  steps.  In  the  first  step  it  is  assumed  that  all 
the  states  of  the  system  are  available  for  feedback  and  a  constant  gain  is  found  such  that 
a  cost  function,  consisting  of  a  weighted  sum  of  the  state  and  control  energy,  is  minimized. 
Often  all  the  states  are  not  available  for  feedback  in  which  case  they  must  be  estimated. 
In  the  second  stage  a  steady  state  Kalman  filter  (or  any  other  observer)  is  designed  to 
estimate  the  states  given  a  set  of  measurements.  The  design  of  the  Kalman  filter  requires 
the  specification  of  process  and  measurement  noise  statistics.  For  the  LQG .  controllers 
used  here  it  was  assumed  that  both  the  process  and  measurement  noise  are  white.  The 
intensity  of  the  process  noise  was  assumed  to  be  the  identity  matrix,  and  the  intensity 
of  the  measurement  noise  was  taken  as  a  constant  multiplied  by  the  identity  matrix.  This 
constant  was  used  to  adjust  the  bandwidth  of  the  Kalman  filter.  It  was  found  experimentally 
that  the  bandwidth  of  the  Kalman  filter  had  to  be  in  the  range  of  1-2  uT  for  satisfactory 
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performance.  For  larger  bandwidths  the  LQG  controllers  were  very  sensitive  to  noise  and 
destabilized  the  compressor.  The  design  of  the  Kalman  filter  also  requires  the  specification 
of  the  cross  correlation  between  the  process  and  measurement  noise,  and  this  was  taken 
as  zero.  This  is  a  crude  assumption  as  there  will  be  correlation  between  the  upstream 
measurement  noise  and  process  noise.  However,  because  the  LQG  controller  was  not  used 
for  its  optimal  properties  this  was  deemed  acceptable. 

It  was  noted  in  Section  3.3.3  that  approximately  4 n  harmonics  were  necessary  to  locate 
the  poles  of  the  first  n  modes.  A  complete  set  of  poles  and  transmission  zeros  of  the 
compression  system  at  <^sta 11  —  1%  is  shown  in  Figure  5.1.  The  poles  and  zeros  corresponding 
to  the  actuator  dynamics  are  not  shown  though  they  are  included  in  the  design  model. 

In  this  model  32  harmonics  were  used  to  compute  the  steady  flow,  and  16  harmonics 
were  included  in  the  state  space  description.  Also  shown  in  the  figure  are  circles  (with 
centers  at  the  origin)  with  radii  of  2  and  3wr.  There  are  five  modes  within  the  2  uiT  circle; 
these  modes  dominate  the  dynamics  at  frequencies  of  practical  interest.  The  other  modes 
are  either  further  into  the  left  half  plain,  that  is,  they  have  fast  dynamics  so  that  they  can 
be  approximated  by  constants  at  low  frequencies,  or  the  modes  closer  to  the  iu— axis  are 
close  to  zeros  and  will  have  little  resonance.  We  also  note  that  there  are  several  poles  and 
zeros  that  almost  cancel.  This  near  pole-zero  cancellation  was  also  observed  at  other  flow 
coefficients  and  distortion  magnitudes. 

With  the  full  set  of  modes  available  the  state  space  model  can  be  reduced  so  that  only 
the  dominant  modes  are  retained.  The  approach  followed  here  was  to  convert  the  state 
space  description  into  modal  form  and  keep  only  the  modes  included  in  the  3  wr  circle, 
typically  15  modes.  It  was  found  that  the  additional  modes  are  necessary  to  capture  the 
dominant  zeros.  If  only  the  first  four  or  five  modes  were  included  the  zeros  were  off  by  as 
much  as  50%.  By  keeping  all  the  modes  inside  the  3  o»r  circles  the  gains  and  phases  were 
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Rea  I  (A)  [u;r] 


Figure  5.1:  Poles  and  transmission  zeros  at  mass  flow  i^aii  —  1%.  The  two 
semicircles  have  radii  2  and  3tvr  respectively. 


approximated  very  well  up  to  2  u>T,  the  bandwidth  of  the  AGV  motors.  This  reduced  order 
model  was  then  used  in  the  LQG  design.  The  modes  of  the  LQG  controller  again  separated 
into  slow  and  fast  modes  so  that  model  reduction  was  also  applied  to  the  controller;  the 
final  controllers  had  between  8  and  10  states.  We  note  that  the  LQG  controller  is  a  dynamic 
controller  while  the  controllers  discussed  below  are  constant  gain  controllers. 

We  note  that  an  LQG  controller  is  a  full  multi-input  multi-output  dynamic  controller, 
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whereas  all  the  remaining  controllers  are  constant  gain  (no  dynamics)  controllers. 


5.2.2  Harmonic  Feedback 

The  second  controller  tested  was  used  by  Haynes  [20]  to  control  rotating  stall  with  uniform 
flow  on  the  same  compressor.  For  uniform  flow  the  modes  have  a  one-to-one  correspondence 
with  the  harmonics  and  go  unstable  one  at  a  time.  This  behavior  was  used  by  Haynes  [20] 
to  experimentally  find  a  set  of  independent  optimal  (in  the  sense  of  achieving  lowest  stalling 
flow  coefficient)  single-input  single-output  controllers  for  the  first  three  harmonics.  Because 
of  the  one-to-one  correspondence  between  modes  and  harmonics  each  controller  acted  on 
only  one  harmonic.  We  will  refer  to  this  as  harmonic  feedback  (HF).  Because  each  har¬ 
monic  feeds  back  only  to  itself  the  (multi-input  multi-output)  HF  controller  has  a  diagonal 
structure.  The  zeroth  harmonic  was  not  considered  important  for  uniform  flow  and  was  not 
used  by  Haynes  [20]. 


5.2.3  Harmonic  Feedback  with  Coupling 

As  shown  in  Chapter  4  distortion  introduces  coupling  between  the  harmonics  and  one  can 
argue  that  the  coupling  can  be  used  to  improve  the  performance  of  a  controller.  For  example, 
the  first  harmonic  can  be  fed  back  to  the  zeroth  and  second  harmonics,  as  well  as  itself.  If 
we  use  only  the  first  three  harmonics  (plus  zeroth),  we  need  to  find  16  gains  and  15  spatial 
phases  (there  is  no  spatial  phase  shift  on  the  zeroth  harmonic).  It  is  very  time  consuming 
to  find  all  these  gains  and  phases  experimentally,  so  it  was  decided  to  employ  coupling 
between  the  zeroth,  first  and  second  harmonics  only.  The  use  of  first  and  second  harmonics 
was  based  on  the  strong  presence  of  these  two  harmonics  in  the  first  and  second  modes 
(see  Chapter  4  for  numbering  of  the  modes).  In  the  experiments  with  the  distortion  it  was 
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found  that  the  compression  system  was  extremely  sensitive  to  changes  in  the  mean  mass 
flow  (that  is,  the  rate  at  which  the  throttle  was  closed,  see  Section  2.4.2)  so  that  it  was 
argued  that  coupling  between  the  zeroth  and  first  harmonics  might  have  a  stabilizing  effect. 


With  this  scheme  the  form  of  the  controller  is 
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kpq  =  \kpq\e'0>\  p,  g  =  0, 1,2, 3. 


The  optimal  values  determined  by  Haynes  were  used  for  fcn,  k22>  and  &33>  and  the  other 
values  ifcoo,  &oi>  &2i  and  were  determined  experimentally.  This  controller  will  be  referred 
to  as  harmonic  feedback  with  coupling  (HFC). 


5.2.4  Distributed  Feedback 

Instead  of  decomposing  the  velocity  perturbations  into  a  sum  of  harmonics  and  feeding 
back  the  harmonics  as  we  did  with  HF  and  HFC,  we  can  view  the  velocity  perturbation  as 
a  spatially  distributed  error  signal  which  is  fed  back  after  it  has  been  multiplied  by  a  gain. 
Also,  because  the  perturbation  is  rotating  around  the  annulus  it  may  be  advantageous  to 
include  a  spatial  phase  shift  0  to  provide  phase  lead  to  compensate  for  the  computational 
delay  and  actuator  dynamics.  These  two  simple  arguments  lead  to  the  control  law 

Sj(e)  =  kS<P(e-0).  (5.1) 

This  feedback  law  will  be  referred  to  as  distributed  feedback  (DF).  The  constant  gain  k  and 
spatial  phase  shift  0  was  determined  experimentally. 
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It  is  interesting  to  compare  this  new  control  law  with  the  HF  controller.  We  can  rewrite 
Equation  (5.1)  in  terms  of  the  individual  harmonics 

8%  =  ke~in08$n  (5.2) 

that  is,  the  nth  harmonic  8<pn  is  rotated  through  an  angle  n/3  and  multiplied  by  a  gain  k 
to  give  nth  harmonic  Sj.  In  contrast,  each  harmonic  of  the  HF  controller  has  its  own  gain 
knn  and  spatial  phase  shift  /3nn 

8%  =  knne-i/3™8$n.  (5.3) 

Another  difference  between  the  two  control  laws  is  that  the  spatial  phase  shifts  of  the  higher 
harmonics  of  the  DF  controller  are  constrained  to  be  multiples  of  that  of  the  first  harmonic. 
As  a  result,  the  phase  relation  that  exists  between  the  harmonics  of  8<f>  is  preserved  by  the 
DF  controller  but  not  by  the  HF  controller. 

5.2.5  Controller  Implementation 

All  the  constant  gain  controllers  can  be  implemented  in  a  computationally  efficient  way 
without  computing  the  harmonics.  Let  ns  be  the  number  of  sensors  around  the  annulus  and 
Fna  be  the  DFT  matrix  corresponding  to  n3  sensors.  Then  we  can  compute  the  harmonics 

84>  =  Fn,8<j>(0 >)  (5-4) 

where  8<f>(9s)  is  a  vector  of  velocity  perturbation  measurements  at  sensor  locations  0\ , . . . ,  81ls 
around  the  annulus.  The  harmonics  of  the  control  signal  are  now  computed  as 

57  =  K84> 

=  KFns8<t>{6a)  (from  5.4)  (5.5) 
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Table  5.1:  Gains  and  spatial  phases  of  the  HF,  HFC,  and  DF  controllers. 
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where  K  is  any  constant  gain  controller.  Let  F~}  be  the  inverse  DFT  matrix  for  synthesizing 
the  AGV  deflections  6j{8&)  at  na  actuator  locations  0i, . . . ,  0n»  around  the  annulus,  then 

*7(*a)  =  F~lS=r 

=  F^KFn,  6<f>(ds)  (from  5.5) 

=  K'5<t>{6s).  (5.6) 

All  the  quantities  are  known  so  the  matrix  K'  can  be  precomputed  and  therefore  we  can 
compute  the  individual  actuator  signals  directly  from  the  individual  sensor  signals  without 
computing  the  harmonics.  If  the  number  of  sensors  is  different  from  the  number  of  actuators 
or  if  the  sensors  and  actuators  axe  at  different  locations  around  the  annulus,  K'  is  not  a 
diagonal  matrix  even  if  the  controller  K  is  diagonal.  The  nondiagonal  form  is  due  to  the 
interpolation  of  the  inverse  DFT. 


5.3  Active  Control  Experimental  Results 

5.3.1  LQG  Controller 

The  stalling  flow  coefficient  for  the  LQG  controllers  are  0.455  and  0.475  for  the  0.8  and 
1.9  dynamic  head  distortions  for  an  increase  in  stable  operating  range  of  1.5%  and  1.1% 


122 


respectively.  The  stable  operating  range  has  been  increased  by  a  small  amount,  but  the  point 
is  that  controllers  based  completely  on  the  model  can  be  used  to  stabilize  the  compressor. 

Because  LQG  controllers  assume  that  certain  phase  relations  exist  between  the  har¬ 
monics  it  is  interesting  to  determine  how  sensitive  the  stalling  flow  coefficient  is  to  changes 
in  the  phase  relations.  This  was  determined  experimentally  by  turning  the  screen  away 
from  the  position  for  which  the  controller  was  designed.  The  resulting  stalling  flow  coeffi¬ 
cients  are  shown  in  Figure  5.2.  We  see  that  for  both  positive  and  negative  offsets  a  point  is 
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Figure  5.2:  LQG  sensitivity  to  offsets  in  distortion  location.  The  dashed  line 
shows  the  stalling  flow  coefficient  without  control. 


reached  beyond  which  the  stalling  flow  coefficient  increases  suddenly  above  that  of  the  open 
loop  stall  point  so  that  the  LQG  controller  has  a  destabilizing  effect.  This  clearly  shows 
the  sensitivity  of  the  LQG  design  to  changes  in  system  parameters. 

The  sensitivity  to  an  offset  in  the  distortion  screen  can  be  explained  as  follows.  Rotating 
the  screen  away  from  its  nominal  position  is  equivalent  to  rotating  both  the  sensors  and 
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actuators  in  the  opposite  direction  through  the  same  angle.  We  will  use  the  latter  approach. 
If  the  angle  through  which  the  sensor  and  actuators  are  rotated  is  90 ,  that  is,  the  screen  is 
rotated  by  —60,  then  the  rotated  input  and  output  harmonics  are 

Sft=6  %e+in0° 

S  ftn°  =  5$ne+in6° 


respectively.  Thus,  the  transfer  function  corresponding  to  the  rotated  configuration  is  ob¬ 
tained  from  the  nominal  transfer  function  G  (see  Equation  (4.15))  as 


diag[l,  e+l6° 

,  e+i26°,  e 

+i3d°]  G  diag[l,  e~l6°,  . 
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The  transfer  functions  on  the  diagonal  axe  unchanged,  but  the  off-diagonal  (or  coupling) 
transfer  functions  are  phase  shifted.  The  LQG  controller  is  a  multi-input  multi-output 
controller  and  uses  all  the  relations  between  the  different  input  and  output  harmonics, 
including  coupling  between  different  harmonics.  By  rotating  the  screen  we  have  changed 
the  system  for  which  the  controller  was  designed,  and  the  phase  change  from  the  coupling 
was  enough  to  destabilize  the  system.  The  sensitivity  of  the  controller  to  changes  in  the 
coupling  between  the  harmonics  is  not  unique  to  the  LQG  controller.  Any  controller  that 
directly  relies  on  the  coupling  will  be  sensitive  to  offsets  in  the  distortion. 
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5.3.2  Harmonic  Feedback  and  Harmonic  Feedback  with  Coupling 


One  question  posed  at  the  beginning  of  the  chapter  was  whether  HF  controllers  designed 
for  uniform  flow  can  stabilize  the  system  in  the  presence  of  distortion.  The  answer  to  this 
question  is  affirmative.  The  stalling  flow  coefficients  of  the  HF  controller  are  0.452  and 
0.475  for  the  0.8  and  1.9  dynamic  head  distortions  for  2.2%  and  1.1%  increase  in  stable 
operating  range  respectively.  Because  the  HF  controller  does  not  use  any  coupling  between 
the  harmonics  it  is  insensitive  to  changes  in  the  circumferential  location  of  the  distortion. 

With  the  HF  controller  under  uniform  flow  inlet  conditions,  Haynes  [20]  achieved  an 
increase  of  7.8%  the  operating  range. 

Another  question  is  whether  the  flow  coefficient  can  be  lowered  even  further  if  we 
take  the  coupling  into  account.  The  HFC  controller  was  used  to  address  this.  Spatial 
phases  for  the  gains  k\2  and  i  were  varied  between  0  -  360°.  It  was  found  that  ki2  =  0 
and  &21  =  2ell35°  gave  the  best  result,  lowering  the  flow  coefficient  to  0.448  for  the  0.8 
dynamic  head  distortion,  a  3%  increase  in  stable  operating  range.  A  similar  procedure  was 
followed  to  determine  &oi  and  k\ o  but  it  was  found  that  the  addition  of  these  gains  led  to 
no  improvement.  The  HFC  controller  thus  increased  the  stable  operating  range  36%  more 
than  the  HF  controller  which  had  no  coupling  between  the  harmonics.  This  improvement 
stems  from  knowledge  of  the  coupling  between  the  harmonics  and  thus  the  HFC  controller 
will  also  be  sensitive  to  offsets  in  the  circumferential  location  of  the  distortion.  The  HFC 
controller  was  not  used  for  the  1.9  dynamic  head  distortion. 

5.3.3  Distributed  Feedback  Controller 

The  DF  controller  was  experimentally  tested  under  both  distortions  and  uniform  flow.  In 
all  cases  it  was  found  that  a  spatial  phase  shift  in  the  range  30  <  (3  <  70°  gave  the  best  per- 
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formance.  Unless  stated  otherwise  0  =  30°  was  used.  The  stalling  flow  coefficients  obtained 
with  this  controller  were  0.445  and  0.472  for  the  0.8  and  1.9  dynamic  head  distortions  re¬ 
spectively.  The  corresponding  improvement  in  stable  operating  ranges  are  3.7%  and  1.7%. 
Both  these  flow  coefficients  are  lower  than  the  corresponding  HF  and  HFC  stalling  flow 
coefficients.  The  DF  controller  does  not  introduce  coupling  between  the  harmonics  so  that 
its  performance  is  invariant  under  changes  in  the  circumferential  location  of  the  distortion. 
An  analysis  comparing  the  HF  and  DF  controllers  will  be  given  after  we  have  discussed  the 
effectiveness  of  the  DF  controller.  Finally,  we  note  that  for  all  the  controllers  the  model 
predicts  that  the  system  should  be  stable  at  the  flow  coefficients  where  we  lose  stability. 
Possible  reasons  for  the  loss  of  stability  are  discussed  in  Section  5.4. 

Figure  5.3  shows  the  PSDs  measured  under  closed  loop  operation  using  the  DF  con¬ 
troller.  For  reference  the  open  loop  PSDs  (Chapter  4)  are  also  shown  in  the  graphs.  The 
figures  show  frequencies  only  up  to  1.1  wr  because  most  of  the  energy  in  the  perturbations  is 
below  1  uT.  The  effect  of  control  is  seen  from  the  attenuation  of  the  peaks  of  the  first  mode 
at  0.43  uT  in  the  zeroth,  first  and  second  harmonics  by  15,  23  and  8  dB  respectively.  The 
third  harmonic  did  not  show  much  open  loop  resonance  and  control  appears  to  have  little 
effect.  In  the  PSDs  we  also  see  that  control  had  very  little,  if  any,  effect  on  the  peak  at  1  o>r 
suggesting  that  the  disturbance  corresponding  to  this  frequency  is  (almost)  uncontrollable 
from  the  AGVs.  This  was  also  true  for  the  other  controllers. 

The  magnitudes  of  the  noise  in  the  first  and  second  harmonics  has  increased  significantly 
from  the  corresponding  open  loop  values  for  frequencies  below  0.1  cnr;  the  situation  is  worse 
for  the  second  harmonic  where  the  noise  had  increased  almost  15  dB  for  frequencies  up  to 
0.35  uiT.  The  PSD  is  a  measure  of  the  energy  in  the  signal  at  a  specific  frequency.  In  all  the 
closed  loop  PSDs  we  see  the  largest  values  at  very  low  frequencies  and  at  1  u>T,  indicating 
significant  energy  at  these  frequencies  and  thus  may  be  causes  for  the  system  losing  stability; 
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this  will  be  discussed  further  in  Section  5.4. 


5.3.4  Comparison  between  Harmonic  and  Distributed  Feedback 

We  now  compare  the  HF  and  DF  controllers.  The  PSDs  prior  to  stall  with  the  HF  controller 
look  very  much  like  the  corresponding  DF  PSDs  and  are  not  shown.  To  highlight  the  differ¬ 
ences  between  the  two  controllers  we  show  in  Figure  5.4  the  differences  in  the  magnitudes 
of  the  PSDs,  that  is, 

A|PSDn|  =  |PSD“f|  —  |PSD°f|  for  n  =  0,l,2,3.  (5.9) 

In  these  figures  positive  values  indicate  that  the  DF  controller  attenuated  perturbations 
more  than  the  HF  controller.  For  all  the  harmonics  the  HF  controller  attenuated  per¬ 
turbations  more  than  the  DF  controller  over  the  range  0.4  -  0.6 u;r  (approximately).  For 
the  zeroth,  second,  and  third  harmonics  and  for  frequencies  below  0.3  u/r  the  DF  controller 
attenuated  disturbances  more  than  the  HF  controller. 

Earlier  in  this  section  we  mentioned  that  there  is  a  large  amount  of  energy  in  the 
noise  at  frequencies  below  0.1  uiT  where  the  DF  controller  generally  had  better  disturbance 
rejection  than  the  HF  controller.  It  may  thus  be  that  noise  at  low  frequencies  contributes 
to  the  loss  of  stability  at  higher  flow  coefficients  for  the  HF  controller.  For  the  zeroth, 
first,  and  second  harmonics  the  DF  controller  attenuated  perturbations  more  than  the  HF 
controller  in  the  range  0.6  -  0.8 cur. 

The  differences  between  the  controllers  can  be  explained  by  looking  at  the  predicted 
root  loci  of  the  two  controllers  systems  shown  in  Figure  5.5.  As  the  gains  are  increased 
from  zero  (no  control)  to  100%  (experimental  values)  the  poles  trace  out  the  paths  shown 
by  the  dots.  In  this  figure  we  see  that  the  HF  controller  stabilizes  the  pole  associated  with 
the  first  mode  slightly  more  than  the  DF  controller,  consistent  with  the  greater  attenuation 
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A|PSD2|  A|PSD0| 


Frequency  [wr]  Frequency  [wr] 

Figure  5.4:  Difference  between  PSD  magnitudes  for  HF  and  DF  controllers. 
A|PSD„|  =  ]PSD^f|  -  |PSD£f|. 
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Figure  5.5:  Root  loci  for  HF  and  DF  controllers.  In  this  figure  the  dots  trace  the 
root  locus  as  the  gains  are  varied  from  0%  (no  control)  to  100%  (experimental 
values).  *  =DF  controller,  x=HF  controller. 
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of  perturbations  by  the  HF  controller  in  the  frequency  range  0.4  -  0.6  u)T  seen  in  Figure  5.4. 
The  gain  of  the  HF  controller  on  the  first  harmonic,  which  is  dominant  in  the  first  mode, 
is  66%  larger  than  the  corresponding  DF  controller  gain  while  the  spatial  phase  shifts  are 
almost  equal,  30°  versus  36°. 

The  biggest  difference  between  the  two  controllers  is  seen  in  the  second  mode.  The 
DF  controller  has  moved  the  pole  associated  with  the  second  mode  essentially  horizontally 
into  the  left  hand  plane  while  the  HF  controller  moved  it  almost  vertically,  offering  little 
stabilization.  The  pole  associated  with  the  second  mode  of  the  DF  closed  loop  system  at 
approximately  u>2  —  0.75  is  in  the  range  0.6  —  0.8  u>T  over  which  the  DF  controller  had 
more  attenuation  (Figure  5.4).  This  is  surprising  as  the  gain  of  the  HF  controller  on  the 
second  harmonic,  which  is  dominant  in  the  second  mode,  is  37%  larger  than  that  of  the 
DF  controller.  However,  the  phases  of  the  two  controllers  are  different,  being  91°  and  60° 
respectively.  Thus,  we  conclude  that  the  better  stabilization  of  the  second  mode  by  the  DF 
controller  is  due  to  the  difference  in  the  spatial  phase  shifts. 

The  third  mode  is  also  stabilized  more  by  the  DF  controller.  In  this  case  the  HF  gain  is 
approximately  30%  smaller  than  that  of  the  DF  controller  and  the  spatial  phases,  110°  for 
HF  and  90°  for  DF,  are  different.  The  HF  controller  moved  the  pole  vertically  so  that  we 
again  attribute  the  better  stabilization  to  the  difference  in  the  phases  of  the  two  controllers. 
We  have  seen  in  the  PSDs  that  the  third  harmonic  did  not  show  strong  resonance.  This 
is  in  part  due  to  the  fact  that  the  third  harmonic  is  not  as  strong  in  the  first  mode  as  the 
zeroth,  first,  and  second  harmonics,  and  in  part  due  to  the  exponential  decay  of  the  higher 
harmonics  upstream  of  the  compressor  which  strongly  attenuates  the  third  harmonic.  It  is 
thus  hard  to  determine  the  importance  of  the  third  harmonic. 

In  Section  3.6  we  saw  that  the  zeroth  harmonic  has  a  strong  presence  in  the  zeroth 
mode.  It  is  thus  surprising  to  see  that  the  HF  controller  moved  this  mode  further  to 
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the  left  than  the  DF  controller,  even  though  the  gain  of  the  HF  controller  on  the  zeroth 
harmonic  is  zero,  while  that  of  the  DF  controller  is  k  =  3.  The  first  harmonic,  however,  has 
strong  coupling  with  the  zeroth  harmonic.  The  larger  gain  on  the  first  harmonic  of  the  HF 
controller  must  have  affected  stabilization  through  this  coupling.  Thus,  the  model  predicts, 
and  experiments  showed,  that  the  DF  controller  stabilizes  the  compressor  better  than  the 
HF  controller. 

The  model  predicts  that  the  DF  controller  also  performs  better  on  distortions  with 
different  magnitudes  and  extents.  The  results  obtained  so  fax  with  the  DF  controller  are 
promising,  but  additional  experimental  results  are  needed  before  its  superiority  can  be 
established.  However,  the  simplicity  and  ease  of  tuning  make  it  a  very  attractive  alternative 
to  the  HF  and  HFC  controllers. 


5.3.5  Other  Experiments  with  Distributed  Feedback  Controller 

We  end  this  section  by  mentioning  some  variations  of  the  DF  controller  that  were  used 
in  experiments.  In  the  previous  discussion  emphasis  was  put  on  the  spatial  phase  shift 
P  with  little  reference  to  the  gain  k.  The  experiments  showed  there  was  no  measurable 
difference  in  the  stalling  flow  coefficients  for  2.5  <  k  <  4.  For  k  <  2.5  the  stalling  flow 
coefficient  increased  and  for  k  >  4.5  the  AGVs  were  continually  driven  into  saturation 
with  no  decrease  in  stalling  flow  coefficient.  The  value  k  =  3  was  thus  used  for  all  the 
experiments. 

It  was  mentioned  earlier  that  the  compressor  was  very  sensitive  to  changes  in  the  mean 
mass  flow  (obtained  through  changes  in  the  throttle  setting).  It  was  believed  it  might  be  due 
to  the  direct  relation  between  mean  mass  flow  and  the  zeroth  harmonic  and  experiments 
were  conducted  where  the  gain  on  the  zeroth  harmonic  was  varied  while  gains  on  the 
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other  harmonics  were  kept  constant.  No  measurable  changes  were  detected  in  stalling  flow 
coefficients  or  PSDs.  The  root  loci  also  show  little  variation  when  the  gain  on  the  zeroth 
harmonic  is  changed.  The  sensitivity  due  to  changes  in  the  mean  mass  flow  thus  remains 
unresolved. 

The  experimental  rig  had  16  hot-wires  so  that  it  was  possible  to  measure  the  first 
seven  harmonics.  It  was  found  that  feeding  back  more  than  the  first  three  harmonics  (in 
the  DF  controller)  led  to  no  decrease  in  stalling  flow  coefficient.  Due  to  the  exponential 
decay  of  the  upstream  flow  field  the  fourth  harmonic  is  attenuated  more  than  20  dB,  and 
the  fifth,  sixth,  and  seventh  harmonics  even  more,  so  that  the  SNR  of  these  harmonics  will 
be  much  smaller  than  one.  In  addition,  the  model  predicts  that  these  harmonics  do  not 
have  a  strong  presence  in  the  first  and  second  modes  so  that  it  is  difficult  to  determine  if 
they  are  important  for  active  control. 

Experiments  were  also  conducted  with  eight  hot-wires  in  the  usual  upstream  locations 
and  eight  hot-wires  installed  downstream  between  the  first  rotor  and  stator.  The  hot-wires 
were  at  the  same  circumferential  locations  (see  Figure  2.4,  p.41).  In  this  sensor  configuration 
the  DF  controller  was  used  to  stabilize  the  system  using  either  eight  upstream  or  eight 
downstream  hot-wires  and  the  HF  controller  was  used  only  with  the  eight  upstream  hot¬ 
wires.  The  stalling  flow  coefficients  are  listed  in  Table  5.2  for  all  the  controllers. 

For  the  DF  controller,  the  stalling  flow  coefficient  for  using  eight  upstream  hot-wires  is 
higher  than  with  16  hot-wires,  giving  an  increase  of  stable  operating  range  of  1.7%  versus 
3.7%  for  16  hot-wires.  The  degradation  in  performance  is  believed  to  be  a  result  of  the 
coarser  localization  of  the  wave  with  the  smaller  number  of  hot-wires.  The  HF  controller 
performed  poorly  in  this  experiment  and  gave  very  little  stabilization  with  only  eight  hot¬ 
wires.  With  the  eight  hot-wires  behind  the  first  rotor  the  DF  controller  had  a  lower  stalling 
flow  coefficient  than  the  corresponding  upstream  case,  0.450  (2.6%)  versus  0.454  (1.7%). 
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Table  5.2:  Stalling  flow  coefficients.  The  first  line  shows  the  stalling  flow 
coefficients  and  percentage  increase  in  flow  coefficient  without  control.  The 
remaining  lines  list  the  stalling  flow  coefficients  and  percentage  increase  in 
stable  operating  range.  The  HF-8  and  DF-8  indicate  that  only  8  hot-wires 
were  used  with  these  controllers.  Percentages  are  computed  with  respect  to  the 
uniform  flow  stalling  flow  coefficient  without  control,  <f>  =  0.460. 


Controller 

Upstream 

Upstream 

Behind  rotor  1 

1.9  dynamic  head 

0.8  dynamic  head 

0.8  dynamic  head 

No  control 

0.480 

4.3% 

0.462 

0.4% 

0.462  0.4% 

LQG 

0.475 

1.1% 

0.455 

1.5% 

HF 

0.475 

1.1% 

0.452 

2.2% 

HFC 

0.448 

3.0% 

DF 

0.472 

1.7% 

0.445 

3.7% 

HF-8 

0.460 

0.4% 

DF-8 

0.454 

1.7% 

0.450  2.6% 

This  result  was  unexpected  as  the  hot-wires  are  in  the  rotor  wakes  so  the  SNR  is  decreased. 
However,  this  degradation  is  partly  offset  because  the  sensing  is  closer  to  the  source  of  the 
perturbations  so  that  the  signal,  especially  the  higher  harmonics,  will  be  captured  better. 
In  addition,  analysis  of  the  transmission  zeros  showed  that  moving  downstream  increased 
the  frequency  of  the  nonminimum  phase  zero  by  50%.  Thus,  we  should  expect  better 
stabilization. 


5.4  Investigating  Loss  of  Stability 

In  the  previous  section  we  mentioned  that  we  lose  stability  at  flow  coefficients  for  which 
the  extended  model  predicts  the  system  should  be  stable.  In  this  section  we  discuss  experi¬ 
ments  that  were  also  performed  to  asses  some  conjectures  concerning  loss  of  stability.  Two 
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simulations  were  also  performed  to  determine  if  actuator  saturation  and  nonlinear  effects 
were  causing  the  loss  of  stability. 

Gysling  [19]  has  found  that  acoustic  modes  in  compression  systems  can  have  a  desta¬ 
bilizing  effect.  To  determine  if  any  acoustic  modes  exist  in  the  experimental  compressor, 
high  response  pressure  probes  were  installed  in  the  plenum  and  the  compressor  was  run  at 
1500,  1900,  2100  and  2400  revolutions  per  minute.  In  all  the  cases  PSD  analysis  showed 
no  new  peaks  and  the  peak  of  the  first  mode  scaled  with  the  mean  rotor  frequency.  The 
experiments  were  repeated  with  the  exhaust  fan  (see  Figure  2.1  on  p.36)  turned  on;  again 
no  differences  were  detected.  It  thus  appears  that  acoustic  modes  are  not  responsible  for 
the  loss  of  stability. 

Day  [7]  has  found  that  some  compressors  do  not  show  modal  perturbations  prior  to 
stall.  Instead,  the  stall  cell  developed  out  of  a  short  length  scale  perturbation.  To  determine 
if  short  length  scale  perturbations  were  present  eight  hot-wires  were  removed  from  upstream 
and  installed  between  the  first  rotor  and  stator.  The  limited  experiments  conducted  with 
this  sensor  arrangement  showed  no  evidence  of  short  length  scale  perturbations.  However,  it 
was  found  that  the  formation  of  the  stall  cell  was  visible  approximately  one  rotor  revolution 
earlier  by  the  hot-wires  behind  the  rotor. 

In  the  PSDs  of  Chapter  4  we  saw  that  there  is  nonlinear  interaction  between  the 
disturbance  at  1  ojt  and  the  first  mode,  resulting  in  peaks  at  uT  ±  wi.  A  similar  type  of 
nonlinear  interaction  was  also  observed  when  the  system  was  excited  with  a  sinusoidal 
signal  of  frequency  we  by  the  AGVs  —  peaks  appeared  in  the  PSD  at  exactly  ujt  ±  we.  In 
this  experiment  the  level  of  excitation  was  5°,  which  is  well  within  the  linear  operating  range 
of  the  AGVs.  The  linearity  of  the  AGVs  was  verified  by  taking  a  full  set  of  transfer  functions 
at  the  same  mean  mass  flow  with  5°  and  10°  excitation  levels.  The  transfer  functions  were 
in  excellent  agreement  so  we  conclude  that  the  nonlinear  coupling  is  not  due  to  nonlinearity 
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in  the  AGVs  and  must  be  due  to  another  phenomenon. 

Prior  to  the  installation  of  the  distortion  screens  uniform  flow  transfer  functions  were 
measured  with  uniform  flow.  The  model  predicts  that  the  harmonics  are  not  coupled  and 
that  the  zeroth  harmonic  is  not  important  for  uniform  flow,  but  the  full  multi-input  multi¬ 
output  transfer  functions  were  measured.  The  uniform  flow  transfer  functions  are  shown  in 
Figure  5.6.  The  coupling  between  the  different  harmonics,  including  the  zeroth,  is  clearly 
visible.  Note  that  the  peaks  in  the  off-diagonal  transfer  functions  are  at  the  same  frequency 
as  that  of  the  first  mode.  This  behavior  is  similar  to  that  we  have  seen  when  distortion 
is  present,  and  therefore  must  be  due  to  nonlinear  effects.  These  transfer  functions  were 
measured  at  <j)  =  0.460,  the  open  loop  uniform  flow  stall  point.  Uniform  flow  transfer 
functions  measured  at  (f>  =  0.470  showed  no  coupling  between  the  harmonics.  At  this 
higher  flow  coefficient  the  resonant  peak  of  the  first  mode  is  significantly  smaller  than  at 
(f)  =  0.460  and  the  coupling  must  therefore  be  a  result  of  the  larger  perturbation  magnitudes. 
This  coupling  was  not  expected  and  led  to  the  determination  of  the  magnitudes  of  the 
perturbations  necessary  for  nonlinear  effects  to  become  important.  This  is  discussed  in  the 
next  section. 


5.4.1  Nonlinear  Simulations 

Analysis  of  the  actuator  deflections  showed  that  the  AGVs  were  often  driven  into  saturation 
prior  to  stall.  To  determine  if  this  was  causing  the  loss  of  stability  the  closed  loop  was 
simulated  with  the  magnitudes  of  the  actuator  deflections  limited  to  15°,  the  maximum 
value  allowed  in  the  experiments.  These  simulations  showed  that  the  closed  loop  system 
should  be  stable  at  the  flow  coefficients  where  the  experimental  system  lost  stability.  It 
was  also  found  experimentally  that  a  decrease  in  the  maximum  AGV  deflections  up  to 
approximately  20%  has  little  effect  on  the  stalling  flow  coefficient.  Therefore,  the  saturating 
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AG  Vs  is  not  the  limiting  factor. 


In  Chapter  4  we  saw  that  the  envelopes  of  the  measured  velocity  perturbations  were 
not  symmetrically  distributed  about  zero,  suggesting  nonlinear  behavior.  The  PSDs  in  that 
chapter  also  showed  effects  characteristic  of  nonlinear  behavior.  Earlier  in  this  chapter  it 
was  mentioned  that  for  small  levels  of  AGV  forcing  nonlinear  behavior  was  also  observed. 
Nonlinear  coupling  was  also  observed  in  the  uniform  flow  transfer  functions  in  Figure  5.6. 
These  observations  suggest  that  the  magnitudes  of  the  velocity  perturbations  may  be  large 
enough  so  that  terms  ignored  during  linearization  may  be  significant.  In  this  section  we 
determine  through  a  simple  simulation  if  this  is  the  case.  The  analysis  presented  here  is  far 
from  complete  and  was  performed  to  see  if  additional  research  is  necessary  in  this  direction. 
Detailed  nonlinear  analysis  (using  Lyapunov  stability  theory)  was  done  by  Mansoux  [31] 
for  uniform  flow,  but  no  such  research  has  been  done  for  distorted  flow. 

To  estimate  the  effect  of  nonlinearity  terms  up  to  second  order  were  kept  in  the  Tay¬ 
lor  series  expansion  of  Equation  (3.28).  The  expression  for  the  pressure  rise  across  the 
compressor  is 

4(rf>  +  50W™)2  -  Mc4's>.  (5.10) 

We  will  only  analyze  the  homogeneous  system  so  the  AGV  forcing  term  is  not  included  here. 
All  the  other  pressure  balance  equations  were  kept  the  same.  Calculations  were  carried  out 
for  the  0.8  dynamic  head  distortion  at  a  flow  coefficient  1%  above  the  open  loop  stalling 
flow  coefficient. 

It  was  found  that  the  magnitude  | <5<£(r  =  0)|  of  the  initial  velocity  perturbation  above 
which  the  system  goes  unstable  was  \6<f>\  =  0.029.  In  Chapter  4  we  saw  that  the  magnitude  of 
the  pre-stall  perturbations  exceeded  this  value.  Although  this  is  a  highly  simplified  analysis 
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it  makes  plausible  the  case  that  nonlinear  effects  may  be  important  in  determining  loss  of 
stability,  but  more  research  is  necessary  in  this  area. 


5.5  Summary 

In  this  chapter  we  presented  the  first  experimental  evidence  that  active  control  can  be 
used  to  increase  the  stable  operating  range  of  a  compressor  operating  with  circumferential 
inlet  distortion.  The  stable  operating  range  was  increased  by  3.7%  and  1.7%  for  the  two 
120°  distortions  of  magnitude  0.8  and  1.9  dynamic  head  respectively  by  the  distributed 
feedback  controller.  It  was  also  shown  experimentally  that  the  harmonic  feedback  controllers 
optimized  for  uniform  flow  were  effective  in  the  presence  of  these  distortions,  although 
the  corresponding  stall  flow  coefficients  were  higher  than  that  of  the  distributed  feedback 
controller.  By  using  the  coupling  between  the  harmonics  it  was  possible  to  lower  stalling  flow 
coefficient.  LQG  controllers  based  completely  on  the  model  also  stabilized  the  compression 
system  at  open  loop  unstable  points. 

The  reason  for  loss  of  stability  is  not  known,  but  experimental  evidence  and  analysis 
of  the  model  suggest  that  the  pre-stall  perturbations  are  large  enough  for  nonlinear  effects 
to  become  significant. 
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6  Summary 


Analytic  and  experimental  studies  were  carried  out  on  a  low  speed  three-stage  axial  com¬ 
pressor  to  assess  the  use  of  active  control  for  increasing  the  stable  operating  range  of  a 
compressor  in  the  presence  of  circumferential  inlet  total  pressure  distortion.  A  set  of  mov¬ 
able  guide  vanes  was  used  for  actuation.  Steady  state  and  dynamic  response  behavior, 
under  open  and  closed  loop  conditions,  were  examined. 

Linear  quadratic  Gaussian  controllers  based  completely  on  the  model  stabilized  the 
system  at  open  loop  unstable  operating  points.  This  is  the  first  experimental  assessment  of 
the  application  of  modern  control  design  methodologies  to  the  control  of  rotating  stall  in 
the  presence  of  total  pressure  inlet  distortion.  The  range  extension  was  1.5%  and  1.1%  for 
0.8  and  1.9  dynamic  head  distortions  respectively. 

Constant  gain  harmonic  controllers  optimized  for  uniform  flow  were  also  able  to  stabilize 
the  compressor  in  the  presence  of  distortion,  extending  the  operating  range  by  2.2%  and 
1.1%  for  distortions  of  0.8  and  1.9  dynamic  heads  120°  extent  respectively.  Taking  into 
account  the  coupling  between  harmonics  of  the  perturbations  increased  the  stable  operating 
range  by  3%  for  the  0.8  dynamic  head  distortion. 

A  new  control  law,  called  distributed  feedback,  was  introduced  and  was  found  experi¬ 
mentally  to  have  superior  performance  in  the  presence  of  distortion.  This  control  law  has  a 
single  gain  and  spatial  phase  shift.  The  distributed  feedback  controller  increased  the  stable 
operating  range  by  3.7%  and  1.7%  for  0.8  and  1.9  dynamic  head  120°  distortions  respec¬ 
tively.  Analysis  has  shown  that  the  improved  performance  results  from  better  stabilization 
of  the  second  mode  and  better  suppression  of  low  frequency  noise.  The  extended  Hynes- 
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Table  6.1:  Comparison  between  Uniform  and  Distorted  Flow 


Uniform  Flow  Distorted  Flow 

Harmonics  are  decoupled.  Section  3.6  Strong  coupling  between  harmonics,  includ¬ 

ing  zeroth. 

Mode  shapes  are  sinusoidal  and  time  invari-  Mode  shapes  are  nonsinusoidal,  change 
ant.  Section  3.6  with  time  as  wave  propagates  around  the 

annulus. 


Envelope  of  velocity  perturbations  uniform 
around  annulus.  Section  4.2.1 

Disturbance  at  1  uT  present  only  in  first  har¬ 
monics.  Section  4.2.2 

SNR  =  1.25  for  the  first  harmonic.  Sec¬ 
tion  4.2.2 

Stalling  flow  coefficient  =  0.460.  Sec¬ 
tions  4.1.2  and  5.3. 

System  becomes  nonminimum  phase  at  4>  = 
0.464(+0.9%).  Section  3.7. 

7.8%  increase  in  operating  range  with  active 
control.  Section  5.3 


Envelope  of  velocity  perturbations  non- 
uniform  around  annulus. 

Disturbance  at  1  uT  stronger  and  present  in 
all  harmonics. 

SNR  =  0.34,  0.63,  0.29  for  zeroth,  first,  and 
second  harmonics. 

Stalling  flow  coefficients:  0.462  (+0.4%) 
and  0.480  (+4.3%)  for  0.8  and  1.9  dynamic 
head  distortions  respectively. 

System  becomes  nonminimum  phase  at  (f>  = 
0.467(+1.5%)  and  <f>  =  0.486(+5.6%)  for  0.8 
and  1.9  dynamic  head  distortions. 

3.7%  and  1.7%  increase  in  operating  range 
for  0.8  and  1.9  dynamic  head  distortions 
with  active  control. 


Greitzer  model  for  compressor  response  to  inlet  distortion  predicts  that  this  controller  would 
also  perform  better  for  other  distortions  with  different  magnitudes  and  extents. 

Transfer  functions  from  the  actuator  guide  vanes  to  the  hot-wires  measured  at  several 
mass  flows  were  in  good  agreement  with  those  predicted  by  the  extensions  of  the  Hynes- 
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Greitzer  model.  This  is  the  first  experimental  evidence  that  quantitatively  establishes  the 
ability  of  the  model  to  predict  the  dynamics  of  small  perturbations  in  the  presence  of 
distortion.  The  transfer  functions  clearly  showed  the  coupling  between  harmonics  which  is 
characteristic  of  the  modes  that  occur  in  the  presence  of  inlet  distortion.  The  open  loop 
stalling  flow  coefficient  for  the  0.8  dynamic  head  distortion  was  correctly  predicted  by  the 
model,  but  was  over  estimated  for  the  1.9  dynamic  head  distortion.  This  discrepancy  is 
probably  due  to  the  fact  that  the  compressor  pressure  rise  characteristic  is  not  known  over 
the  complete  range  of  the  velocities  present  for  the  large  distortion. 

The  experimental  results  revealed  nonlinear  dynamic  behavior  connected  with  stall 
inception.  The  envelope  of  the  velocity  perturbations  follows  the  trend  predicted  by  the 
model,  but  is  not  symmetrically  distributed  about  the  mean.  Power  density  spectra  revealed 
an  interaction  between  the  first  mode  and  the  disturbance  at  one  rotor  revolution,  and  a 
nonlinear  simulation  showed  that,  at  the  magnitudes  of  perturbations  observed  prior  to 
stall,  nonlinear  effects  may  be  significant. 


Future  Research 

Several  questions  remain  unsolved  and  need  additional  investigation.  The  most  important, 
loss  of  stability  under  active  control,  is  not  understood.  Spectral  analysis  has  shown  that 
all  the  controllers  suppressed  perturbations  of  the  most  unstable  mode  well,  but  there  is 
a  significant  amount  of  energy  in  low  frequency  noise  as  well  as  at  one  rotor  revolution. 
The  noise  far  upstream  and  at  the  throttle  exit  does  not  follow  the  trends  observed  in 
the  experiments,  and  it  is  not  known  if  the  low  frequency  noise  contributes  to  the  loss  of 
stability.  In  addition,  experimental  data  showed  various  forms  of  nonlinear  behavior,  and 
it  is  not  known  whether  it  is  driving  the  system  unstable.  The  disturbance  at  one  rotor 
revolution  is  significantly  stronger  in  the  presence  of  distortion  and  interacted  in  a  nonlinear 
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way  with  the  first  mode.  The  cause  of  this  disturbance  is  not  known,  and  it  is  also  not 
known  if  it  contributes  to  the  loss  of  stability. 

The  extent  and  intensity  of  the  distortions  have  been  assumed  known  and  time  invariant 
but  this  is  not  a  realistic  assumption  under  all  circumstances.  The  effectiveness  of  the 
controllers  must  be  verified  in  the  presence  of  unknown,  time  varying  distortions.  If  the 
rate  at  which  the  distortion  changes  is  slow  compared  to  the  dynamics  of  rotating  stall, 
the  distortion  can  be  approximated  as  steady,  and  it  may  be  possible  to  find  controllers, 
optimized  for  the  instantaneous  distortion,  in  real  time  using  recursive  algorithms. 

Schulmeyer  [43]  showed  that  restaggering  inlet  guide  vanes  by  10°  decreased  the  nonuni¬ 
formity  in  the  flow  by  50%  and  gave  an  increase  in  pressure  rise  of  5.3%.  By  combining 
restaggering  inlet  guide  vanes  with  active  control  we  should  be  able  to  achieve  a  larger 
increase  in  the  stable  operating  range.  If  the  distortion  is  changing  slowly  the  restaggering 
can  be  done  dynamically  in  real  time.  This  combined  approach  is  promising  and  warrants 
further  investigation. 
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Appendix  A  Fourier  Convolution  Matrix 


Let  two  Fourier  series  be  given  by 

nf 

/w =/»+£(/„  cos  p9  +  fsp  sin  p9)  (A.  1) 

p= i 

n9 

g(9)  =9o+  cos  q9  +  gaq  sin  qd).  (A.2) 

<7=1 

Multiplying  the  Fourier  series  f(6)  with  rij  harmonics  with  g{6)  with  ng  harmonics 
gives  another  Fourier  series  p(6)  with  rif  +  ng  harmonics.  Note  that  a  Fourier  series  with  n 
harmonics  has  2n  +  1  Fourier  coefficients. 

We  would  like  to  compute  the  product  p(9)  =  f{9)g{9)  and  write  it  in  the  form  p  = 
F(f)g  where  F(f)  is  the  (2 (n/  +  ng)  +  1)  x  (2 ng  +  1)  Fourier  convolution  matrix  with 
elements  that  are  functions  of  the  coefficients  of  f(9),  g  is  a  2 ng  +  1  vector  with  elements 
the  coefficients  of  g{9)  stored  in  the  format 

9  =  [ffcO>  9cli  9sli  9c2i  5s2>  9cngi  ^sn9]^)  (A. 3) 

and  p  is  a  2(n/  +  ng)  + 1  vector  with  elements  the  coefficients  of  the  product  p{9)  stored  in 
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a  format  similar  to  that  of  g.  Write  the  product  as 

ns 

f(9)g(9)  =  f0go  +  /o  X  cos  $  +  9sg  sin  qO)  (A.4) 

9=1 

nf 

+  Yl  ( fcp  cos p9  +  fap  sin p9)  gp  (A.5) 

p=1 

71  f  rig  n/  Tig 

+ X  /<* cos  p(>  X 5c? cos  761 + X  f°p cos  X  9s<> sm  qd  (A-6) 

p=i  9=1  p=i  9=1 

n/  Tig  Tif  ng 

+  X  A?  Sm?e  X  5c<?  C0S  9°  +  X  Ap  sin^  X  9s<! sin  90  (A-7) 

p- 1  9=1  p=l  9=1 

The  coefficients  of  the  expressions  in  (A.4)  and  (A.5)  are  easy  to  handle.  The  expressions 
in  (A. 6)  and  (A.7)  have  products  of  the  form  cos p9 cos g0,  cos p0sing0,  sinp0cosg0,  and 
sin  p0  sin  q9.  From  standard  trigonometric  identities  these  are  given  by 

fcpgcq  cos p9 cos  q9  =  [cos (p  +  q)9  +  cos (p  -  q)9]gcq  (A.8) 

/cP5s?  cosp0sing0  =  ^[sin(p  +  g)0  -  sin(p  -  q)9]gsq 

=  ^[sin(p  +  ?)0  -  sign(p  -  q)  sin(|p  -  9|)0]pS9  (A.9) 

fspgcq  sin p9  cos  q9  =  ^ [sin (p  +  g)0  +  sin(p  -  g)%C9 

=  [sin(p  +  q)9  +  sign(p  -  g)  sin(|p  -  q\)9]gcq  (A.10) 

fspgsq  sin p9 sin  q9  =  -/sp2[cos(p  +  q)9  -  cos (p  -  g)0]gs<j  (A.ll) 

where  |  •  |  denotes  absolute  value  and 

+1  if  p  >  g, 

sign (p  -  9)  =  <  0  if  p  =  9,  (A-12) 

—1  ifp<g. 
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The  sign  and  absolute  value  modifications  are  included  to  simplify  software  implementation. 
Therefore,  for  each  of  the  four  products  and  for  each  combination  of  p  and  q,  we  get  two 
entries  in  the  matrix:  one  corresponding  to  p  +  q  and  one  to  p  —  q.  This  brute  force  method 
is  not  the  most  efficient  way  to  construct  the  Fourier  convolution  matrix  but  works  quite 
well  for  practical  sized  problems. 

We  illustrate  the  use  of  the  Fourier  convolution  matrix  with  an  example.  Let 

f{9)  =  1  +  2cos0  +  4sin0 

g(9)  =  —  2  +  6cos0  —  8sin0,  (A.13) 


then  p(9)  =  — 12  +  2  cos  9  —  16  sin  9  +  22  cos  20  +  4  sin  29.  Vectorizing  we  get 


Hf)  = 


( ,  !  2' 


2  1  0 

4  0  1 

0  1  -2 


0  2  1 


(A- 14) 


and  g  =  [—2,6,  —  8]r.  Multiplying  the  matrix  F(f)  with  the  vector  g  we  get 


p  =  [—12, 2,  —16, 22,  (A. 15) 


The  elements  ofp  are  exactly  the  coefficients  of  p(9). 

This  example  also  illustrates  the  increase  in  the  number  of  harmonics  that  we  get  when 
we  multiply  two  Fourier  series.  For  software  implementation  we  need  to  take  care  of  this 
increase  in  the  number  of  harmonics.  This  is  especially  true  for  the  least  squares  minimiza¬ 
tion  routine  that  uses  an  iterative  algorithm  to  solve  for  the  steady  flow.  The  approach 
adopted  here  was  to  truncate  the  order  of  the  Fourier  series  after  every  multiplication  to  a 
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predetermined  number  of  harmonics.  The  number  of  harmonics  that  is  needed  was  found 
through  simulations  and  is  discussed  in  Section  3.3.3. 
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Appendix  B  AGV  Dynamics 


Expressions  for  the  pressure  rise  across  and  flow  perturbations  through  nonuniformly  stag¬ 
gered  AG  Vs  for  uniform  steady  flow  were  derived  by  Longley  [26].  In  this  appendix  we 
extend  the  results  to  the  nonuniform  steady  flow  case.  In  contrast  to  the  chapters  in  the 
thesis  the  derivation  in  this  appendix  uses  absolute  values  of  some  variables  and  we  will 
nondimensionalize  the  expression  as  we  proceed.  In  particular,  the  variables  6a,  pt,  and  x 
axe  not  nondimensionalized;  the  remaining  variables  have  their  usual  meaning. 

Consider  two  vanes  at  circumferential  locations  9  and  9  +  A 9  with  deflections  7  and 
7  +  A7  respectively  as  shown  in  the  figure  below.  The  widths  between  the  vanes  at  station 
2  (inlet)  and  station  3  (outlet)  are  given  by 

e 

e  +  Ae 

e 


wm  =  rA  9  —  ^-[(7  +  A7)  —  7] 
w(3)  =  fA9  +  y  [(7  +  A7)  -  7] 
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respectively.  Dividing  through  by  rA9  and  talcing  the  limit  A9  ->  0  we  get 


wm 

Tab 


i 


w<*> 

Tab 


«  .  ba  i 
=  1  +  TrJ 


where  we  used  a  prime  to  indicate  derivative  with  respect  to  9.  Mass  conservation  requires 
that  wm<f>m  =  w{3)4>(i\  so  that,  to  first  order, 


(j)m  =  (1  +  Ma7')<£(3) 


(B.l) 


where  we  have  used  the  definition  /xa  =  ba/r.  Linearizing  this  equation  about  the  nonuniform 
steady  flow  and  steady  deviation  we  get 

6<p(2)  =  (1  +  /ia7s  ')HW  +  Ha.fcSj'.  (B.2) 


Next  we  find  an  expression  for  the  total  pressure  across  the  AG  Vs.  Along  a  streamline 

(B.3) 


Pt  dip 

- j-  — -  =  constant 

p  dt 


cx  —  <puw 

_  dp 
dx  * 


(B.4) 


Assuming  a  linear  axial  velocity  distribution  between  stations  (2)  and  (3)  we  can  write 


c,(x)  =  C«>  +  f  (4»  -  4”). 
ba. 

With  this  velocity  relation  we  have  from  Equation  (B.4) 

<p{x,  t)  =  c(t)  +  xc(2)  +  |^(43)  -  42)) 

=  J_£c(t)  df»_ 

uw  dt  uw  dt  dt  2 ba  dt 
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Prom  Equation  (B.3)  we  get 


Pt  -Pt  _  dPm  _  d(P(3) 
p  dt  dt 


Simplifying  and  nondimensionalizing  this  equation  gives 


D<3)  _  0<2>  d  1 

Z-P-  =  -  #”))■ 

pu^  OT  2 


Prom  Equation  (B.l)  we  get  the  first  order  approximation 


IfoW  _*(*>)*(!  +  \w')<t>™. 


Substituting  this  into  Equation  (B.5)  and  simplifying  gives 

(3)  _  (2)  g  i 

Linearizing  this  equation  about  the  nonuniform  steady  flow  gives 

5p^  —  Sp[2)  .  1  t\rl(3)  Ma  jl  j :•/ 

— - 2^—  =  “Ma(l  +  xMa7  )H  ~  ^<t> s<*7  • 

put  2  2 


Equations  (B.2)  and  (B.6)  are  the  equations  needed  for  the  linearized  dynamics. 
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Appendix  C  Transfer  Functions 

C.l  Complex  Input-Output  Transfer  Functions 

In  this  appendix  we  derive  expressions  to  compute  the  complex-input  to  complex-output 
harmonic  transfer  functions  given  the  real-input  to  real-output  harmonic  transfer  functions. 
The  complex-input  to  complex-output  transfer  functions  will  be  referred  to  as  complex 
transfer  functions. 

For  simplicity  we  derive  the  results  for  one  harmonic  plus  the  zeroth  harmonic  and 
use  the  subscripts  c,  s,  and  0  to  refer  to  the  cosine,  sine  and  zeroth  harmonic  coefficients 
respectively.  The  results  hold  for  all  harmonics  n  >  0  and  the  input  and  output  harmonic 
numbers  do  not  have  to  be  the  same.  Consider  the  real-input  real-output  transfer  function 
defined  by 


Vo(*w) 

9oo(iu)  goc(iu)  goa(w) 

u0(iu) 

&(«*>) 

= 

5co(M  gcc(iw)  9cs(iu) 

uc[iu)) 

- 1 

3 

• 

£>t 

9so{iu)  gsc{iu)  ffss(fw) 

ua(iu ) 

To  keep  the  notation  uncluttered  we  will  omit  the  dependency  of  the  transfer  functions  on 
iuj  in  most  of  the  expressions.  Define  the  complex  harmonic  output 


y{t)  =  yc(t)  +  iys(t) 

y{iuj)  =  yc(w)  +  iya{iu) 

(C.2) 

y*(t)  =  Vc(t)  -  iys{t ) 

y(iu>)  =  yc{iu)  -  iys{iu) 

(C.3) 

'A* 

where  y*  is  the  complex  conjugate  of  y.  Note  that  y*(iu)  ±  y(w).  Similarly,  for  the  complex 
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harmonic  input  u  we  have 

u(t)  =  uc(t)  +  ius(t)  u(iu)  =  uc(iu)  +  ius(iu ) 

u*(f)  =  uc(t )  —  iua(t)  u(iuj)  =  uc(iu)  —  ius(iu}). 

As  noted  in  Section  C.l  the  definitions  of  the  complex  harmonics  do  not  correspond  to  the 
standard  relation  between  real  and  complex  Fourier  coefficients.  The  two  variables  y*  and 
u*  are  introduced  to  simplify  the  similarity  transformation.  The  similarity  transformation 
matrix  and  its  inverse  that  corresponds  to  these  definitions  are 


T  = 


1  0  0 
0  1  i 
0  1  -i 


■7-1  _ 


2  0  0 
0  1  1 
0  —i  i 


Applying  the  similarity  transformation  to  Equation  (C.l)  we  get  the  complex  transfer  func 
tions 

t  r 

9  oo  \  (9oc  —  J9os)  j(5oc  +  iffos) 

9co  +  l5so  |[(ffcc  +  9ss)  ~  *(Pcs  —  9sc)\  j[(5cc  —  9ss)  +  i(gcs  +  9sc)} 

9 co  ^9so  2  [(^cc  9ss)  ^(*?cs  d"  5sc)]  2^9 cc  9ss)  d*  i{9cs  5sc)] 

(C.4) 


Vo  fa) 

y{icj) 

= 

" 

uQ(iuj) 

u(iu)) 

u(iw) 

9oo(iv)  gou(iu)  gQ*(i u) 
gy0{iui)  gyu(iui) 

9'y  0M  5*U(^) 


u0(iu>) 

u(iu>) 

u(icj) 

(C.5) 


We  are  especially  interested  in  the  transfer  functions  g00(iu>),  gou(iu),  gy0(iu),  and  gyu(iuj). 
For  positive  values  of  u  these  give  the  complex  transfer  functions  between  the  different 
harmonics  and  a  peak  in  the  magnitude  of  the  transfer  function  will  indicate  a  wave  travel¬ 
ling  in  the  same  direction  as  the  rotor.  Because  the  input  and  output  signals  are  complex 
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these  transfer  functions  do  not  satisfy  the  complex  conjugate  symmetry  property  and  the 
transfer  functions  at  negative  frequencies  look  completely  different  from  those  at  positive 
frequencies.  Therefore,  if  we  want  to  have  a  complete  picture  of  the  system,  we  need  to 
look  at  both  positive  and  negative  frequencies.  Fortunately  we  do  not  have  to  re-evaluate 
the  transfer  functions  at  negative  frequencies  (which  is  computationally  expensive)  —  they 
can  be  computed  from  the  real  transfer  functions  by  using  a  few  simple  relations  that  we 
derive  next. 

From  Fourier  transform  properties  we  know  that  the  Fourier  transforms  of  real  val¬ 
ued  functions  have  complex  conjugate  symmetry,  that  is,  if  f(t)  is  a  real  function  with 
Fourier  transform  /(iw)  then  Applying  this  to  the  transfer  function  of 

Equation  (C.l)  we  get 


ffccM  =  0cc(-*w) 


9ls(iuj)  =  0cs(-«*>) 


Using  these  relations  we  find 

0ou(-*w)  =  ^bocM  -  *0osM]  (C-6) 

gy0(-iu)  =  g*0(iu)  +  ig*0(iu)  (C.7) 

gyu{-iu)  =  ^[(SccM  +  gtsiiu))  -  -  ffs*c(w)].  (C.8) 

Equations  (C.4),  (C.5),  and  (C.6)-(C.8)  are  all  we  need  to  compute  the  complex  transfer 
functions  at  all  frequencies  given  the  real  transfer  functions  at  positive  frequencies. 

Other  interesting  relations  can  be  derived  by  applying  the  conjugate  symmetry  to  the 
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other  transfer  functions  in  Equation  (C.4).  The  following  relations  are  easy  to  derive 

gyu(-iuj)  =g*„(icu) 

C.2  Measurement  of  Real  Transfer  Functions 

Various  time  and  frequency  domain  methods  exist  for  identifying  multivariable  systems; 
see  for  example  [45].  A  major  consideration  for  measuring  the  compression  system  transfer 
functions  was  the  quality  of  the  measurements  at  low  signal  to  noise  ratios.  Of  all  the 
available  methods  the  nonparametric  frequency  domain  correlation  method  is  considered  the 
most  robust  method  for  measuring  transfer  functions  in  noisy  environments  (see  Rake  [40] 
and  Wellstead  [47]).  In  addition,  it  is  also  robust  against  even  ordered  nonlinearities  and 
methods  exist  to  eliminate  odd-ordered  nonlinear  effects  [11].  It  is  interesting  to  note  that 
this  is  one  of  the  oldest  identification  methods  and  has  been  available  in  commercial  transfer 
function  analyzers  since  the  Sixties  [9].  We  first  present  the  method  for  a  single-input  single¬ 
output  system  and  then  we  generalize  it  to  the  multivariable  and  closed  loop  cases. 

Consider  a  single-input  single-output  system  with  transfer  function  g(s)  as  shown  in 
Figure  C.l.  If  the  input  to  the  system  is  a  cosine  with  constant  amplitude  uq  and  frequency 
uq  the  output  y(t)  of  the  system  is 

y(t)  =  |p(iu;o)|  cos(u/ot  -1-  <p)uo  +  n(t)  +  transients 

where  n(t)  represents  the  noise  in  the  measurement.  If  the  system  is  stable  the  transients 
will  go  to  zero  and  the  system  will  reach  steady  state  —  we  assume  that  this  is  the  case. 
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Uo  COS  UJO  t 


y 


Figure  C.l:  Open  loop  transfer  function  measurement. 


Multiplying  the  output  by  cos  u>ot  and  integrating  over  k  periods  to  =  2tt/u)q  we  get 


2  r^t  o 

Pr  =  —  /  y(t)  cos(u;0t)dt 
kto  J0 

rkto 

=  \g(iuo)\cos{ip)u0/2  + t—  n(t)  cos  (w0*)dt.  (C.9) 

kto  Jo 


If  the  noise  is  zero  mean  and  independent  from  the  input  cosine  or  if  the  noise  covariance  is 
bounded  by  an  exponential  function  (a  weak  assumption  for  stationary  processes)  the  last 
integral  goes  to  zero  as  k  goes  to  infinity  [45].  Similarly,  forming  the  product  between  the 
output  and  sin  wot  we  get 


1  fkt° 

Pi  =  r~r  /  y(t)  sin(u;0t)dt 
kto  Jo 

1  fkt° 

=  |p(io>o)|  sin(<p)uo/2  +  —  /  n(t)  sin(wot)dt.  (C.10) 

kto  Jo 


From  Equations  (C.9)  and  (C.10)  we  can  compute  the  magnitude  and  phase  of  the  transfer 
function  at  uq  from 


|p(uu0)|2  =  (-)2(p?+p?)  (C.ll) 

no 


<p  =  argp(wo) 


=  arctan(pi/pr).  (C.12) 
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There  axe  various  ways  of  viewing  this  procedure.  One  view  is  that  of  correlating  out 
the  known  signal  at  the  output  of  the  system  as  shown  in  Figure  C.2.  An  alternative  view 


COS  CJot 


sin  uq  t 


Figure  C.2:  Transfer  function  measurement:  correlation. 


is  that  of  filtering  the  output  y(t)  with  two  finite  impulse  response  (FIR)  filters  [9]  as  shown 
in  Figure  C.3.  The  filter  impulse  responses  and  transfer  functions  are 


Figure  C.3:  Transfer  function  measurement:  FIR  filtering. 


^r(^)  —  i 

jfcfecos(wof) 

0  <  t  <  t0 

< — >  hT(iuj)  = 

-1  u  Sm(nkuj/uj0)  __i7rfca,/(Jn 

nk  uq  1  —  (w/wo)2 

0 

< 

otherwise 

h{{t)  =  < 

f 

jsi  sm(uj0t) 

0  <  t  <  t0 

i 

5r 

ii 

-i  sin(-rrku)/uJo)  c_i„ku/uJ0 
irk  1  —  (■ w/wq )2 

0 

otherwise 

The  magnitudes  of  the  transfer  functions  corresponding  to  the  impulse  responses  are 
shown  in  Figure  C.4  for  two  values  of  k.  For  the  figures  the  magnitudes  of  the  transfer 
functions  have  been  multiplied  by  two  to  compensate  for  the  factor  two  in  Equation  (C.ll). 
Note  further  that  the  bandwidths  are  quite  narrow  so  that  only  ±10%  of  the  frequency 
axis  about  the  normalized  center  frequency  w/uiq  is  shown.  Only  a  few  of  the  side  lobes 
are  shown  and  the  dotted  lines  in  the  graphs  indicate  the  envelopes  of  the  side  lobes.  If 
we  define  the  bandwidth  Auj/wo  of  the  filter  as  the  width  between  the  first  zeros  of  the 
side  lobes  then  Aw/wo  =  2 /k.  Nondimensionalized  frequencies  in  the  range  0.1  <  u»o  <  3.0 
rotor  revolutions  correspond  to  absolute  frequencies  in  the  range  4-120  Hz,  so  that  for  an 
integration  time  of  10  seconds,  40  <  k  <  1200  and  0.05  <  Auj/uq  <  0.0017.  The  standard 
deviation  of  the  transfer  function  estimate  G  is  (see  [40]) 


STD(G(iw0))  = 


I  PSDn(u,0) 


u0|<3(iwo)|  y  kto 
where  PSDn(uo)  is  the  value  of  the  noise  power  density  spectrum  at  ujq. 


(C.13) 


The  method  can  be  generalized  to  the  multi-input  multi-output  case  in  several  ways. 
One  way  is  to  excite  one  input  at  a  time  and  measure  all  the  outputs,  that  is,  we  are 
measuring  a  single-input  multi-output  transfer  function.  By  using  linearity  we  can  measure 
the  complete  transfer  function  by  repeating  the  procedure  for  the  other  inputs  one  at  a  time. 
A  significant  amount  of  time  can  be  saved  by  exciting  an  input  with  different  frequencies 
(and  appropriately  selected  phases)  at  the  same  time.  Alternatively,  we  can  excite  the 
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different  inputs  with  different  frequencies.  Procedures  exist  for  constructing  these  signals, 
called  multisines,  so  as  to  minimize  the  maximum  amplitude  of  the  excitation  signal,  see  [13]. 
A  disadvantage  of  multisines  is  that  the  control  power  needed  goes  up  with  the  square  of  the 
number  of  sines.  For  the  current  actuators  it  was  found  that  with  two  or  more  frequencies 
at  a  time  the  AGV  motor  currents  exceeded  their  maximum  values  and  were  electronically 
limited.  This  was  considered  undesirable  as  it  causes  distortion  of  the  excitation  signal  so 
that  only  one  frequency  and  one  input  was  used  at  a  time. 

It  is  possible  to  generate  the  input  signal  so  that  the  resulting  wave  travels  in  any 
direction.  Although  it  is  tempting  to  generate  a  wave  rotating  in  the  same  direction  as  the 
rotor  it  was  found  that  this  would  easily  destabilize  the  system  for  excitation  frequencies 
close  to  that  of  the  first  and  second  modes.  Even  if  the  compressor  stayed  stable  the 
perturbations  were  too  large  to  be  considered  linear  operation  so  that  corotating  excitation 
signals  were  avoided  dining  transfer  function  identification  experiments. 

The  above  method  can  also  be  extended  to  closed  loop  transfer  function  measurements 
provided  that  proper  care  is  taken  to  prevent  the  method  from  giving  biased  results.  The 
bias  arises  because,  under  closed  loop  operation,  the  noise  is  being  fed  back  so  that  the 
noise  at  the  input  of  the  system  is  correlated  with  the  noise  at  the  output  of  the  system. 
Unbiased  measurements  can  be  obtained  if  the  system  is  excited  with  an  external  signal  r 
that  is  uncorrelated  with  the  noise  as  shown  in  Figure  C.5.  In  this  figure  K(s)  represent  the 
controller,  M(s)  the  motors  and  amplifiers,  and  G(s)  the  compressor  and  anti-alias  filters. 
We  are  interested  in  measuring  G(s) .  This  is  done  by  measuring  the  transfer  functions  from 
r  to  S'y  and  r  to  5<j>  at  the  same  time  and  recovering  G(s)  from  these  two  transfer  functions. 

The  detailed  closed  loop  system  is  shown  in  Figure  C.6  where  we  indicated  transfer 
function  measurements  to  the  velocity  as  well  as  to  the  plenum  pressure  perturbations. 
Ap(s)  and  A^,(s)  are  known  anti-alias  filters.  First,  assume  the  loop  is  not  closed  so  that 
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Figure  C.5:  Closed  loop  transfer  function  measurement. 


the  transfer  function  from  r  to  the  AGV  deflections  6*y  is 


■<^5) 
6<f> 


Figure  C.6:  Detail  transfer  function  measurement. 


£7  =  Mr 
=  Gir. 


The  open  loop  transfer  function  from  r  to  the  velocity  perturbations  6<fr  is 


Scf>  -  A^G^Mr 


—  A'f.GtjtGir 

=  G2r , 


(C.14) 


(C.15) 

(C.16) 
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and  the  open  loop  transfer  function  to  the  plenum  pressure  perturbations  is 


p<5)  =  ApGpMr 
=  ApGpG\r 
=  Gzr. 

Prom  Equations  (C.14)-(C.18)  we  get 

g^  =  a;1g2gj1 

Gp  =  A;lGzGil. 

Under  closed  loop  operation  the  transfer  function  from  r  to  S-y  is 

67  =  Mr-  MKA'frGfS'y 
=  {I  +  MKA(j>G<i>)-lMr 
=  G\r, 

and  the  closed  loop  transfer  function  from  r  to  6(f)  is 

6(f)  =  A(pG(f)6'y 
=  A^G^Gir 
=  G2r, 

and  the  closed  loop  transfer  function  from  r  to  p£5)  is 

p'5)  =  ApGp51 

=  ApGpGir 

=  Gzr. 


(C.17) 

(C.18) 

(C.19) 

(C.20) 


(C.21) 


(C.22) 

(C.23) 


(C.24) 

(C.25) 
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From  Equations  (C.21)-(C.25)  we  get 


Gt  =  A?G2G? 

(C.26) 

gp  =  ^ig3g:1. 

(C.27) 

We  note  that  the  transfer  functions  G\,  G2,  and  G3  defined  in  Equations  (C.14)-(C.18) 
are  different  from  the  ones  defined  in  Equations  (C.21)-(C.25)  respectively  even  though 
they  have  been  defined  with  the  same  symbols  —  the  former  ones  are  open  loop  transfer 
functions  while  the  latter  ones  are  closed  loop  transfer  functions.  The  form  of  the  final 
expressions  for  G $  and  Gp  axe  exactly  the  same  for  the  open  and  closed  loop  cases  so  that 
the  same  program  can  be  used  to  compute  the  required  open  loop  transfer  functions. 

We  note  an  interesting  feature  of  this  approach  —  it  does  not  require  the  controller 
to  be  “square,”  that  is,  the  number  of  inputs  to  the  controller  5<p  does  not  have  to  equal 
the  number  of  outputs  from  the  controller  S'y.  However,  we  do  require  the  transfer  function 
from  r  to  67  to  be  square.  Measuring  the  transfer  functions  as  outlined  above  has  an 
additional  advantage  that  we  do  not  need  to  know  the  transfer  functions  of  the  individual 
building  blocks,  including  the  transfer  function  of  the  controller.  Everything  that  is  needed  is 
included  in  the  transfer  function  G\  from  r  to  Sj  that  is  measured.  The  transfer  function  G\ 
can  be  measured  accurately  because  the  AGV  deflections  can  be  measured  accurately  with 
the  shaft  encoders.  This  method  further  does  not  require  the  motor-amplifier  channels  to 
be  identical;  even  though  matching  was  not  required  for  the  transfer  function  measurements 
the  individual  amplifiers  were  all  adjusted  so  that  their  bandwidths  matched  to  better  than 
3%  as  a  mismatch  between  the  channels  will  affect  control. 

In  all  the  transfer  function  measurements  the  system  was  excited  for  5  seconds  before 
data  collection  was  started  to  ensure  that  all  transients  died  out  and  that  the  system  has 
reached  steady  state.  The  integration  time  was  set  to  10  seconds  for  all  frequencies  and 
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excitation  frequencies  were  always  selected  so  that  a  full  number  of  periods  would  fit  into 
the  integration  time.  The  magnitude  of  the  excitation  signal  was  varied  between  2°  for 
frequencies  close  to  those  of  the  first  and  second  modes,  to  5°  for  frequencies  fax  away 
from  the  first  and  second  modes.  This  small  level  of  excitation  ensured  linear  operation  of 
the  AGVs  and  kept  the  amplifiers  from  saturating  at  frequencies  up  to  about  three  times 
the  rotor  frequency.  Furthermore,  frequencies  at  one  half,  one,  and  two  times  the  rotor 
frequencies  were  not  used  during  transfer  function  measurements. 
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Appendix  D  System  Identification 

This  appendix  outlines  a  procedure  to  identify  the  state  space  matrices  from  a  set  of  mea¬ 
sured  transfer  functions. 

For  simplicity  the  presentation  is  given  for  single-input  single-output  (SISO)  system, 
and  extension  to  multi-input  multi-output  (MIMO)  systems  is  briefly  commented  on.  A 
brief  description  of  parameterization  is  also  given. 

Let  g(s)  be  the  measured  transfer  function.  We  would  like  to  find  a  corresponding  state 
space  description.  The  procedure  followed  here  is  first  to  fit  a  transfer  function  of  the  form 
b(s)/a(s)  to  the  measured  transfer  function,  and  then  find  a  state  space  description  for  it, 
keeping  in  mind  that  we  are  interested  in  the  MIMO  case. 

Though  we  use  s  —  iu>  as  frequency  variable,  the  method  is  directly  applicable  to 
sampled  data  systems  by  replacing  s  with  z  =  e*n 


D.l  SISO  System  Identification 


We  fit  a  model  of  the  form  b(s)/a(s)  to  the  measured  transfer  function  g(s)  by  minimizing 
the  nonlinear  least  squares  (NLS)  criterion 


c 


LJ 


M 

a(s) 


~9(s) 


2 


(D.l) 


Levy’s  method  (see  [48])  is  often  used  to  turn  the  NLS  problem  into  a  linear  least 
squares  problem  with  cost  function 

c'  =  lb(s)  _  a(s)s(s) I2  (D-2) 

LJ 
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that  is  easy  to  solve.  This  method  often  gives  poor  results  as  we  are  ignoring  a  factor  l/a(s) 
in  the  last  equation;  the  correct  expression  for  the  cost  function  is 

£  -«(*)«(*)]  •  (D.3) 

LJ  '  ' 

Levy’s  method  can  be  improved  by  iterating  the  procedure  as  follows.  We  start  the 
procedure  by  assuming  an  initial  estimate  ao(s)  =  1  and  solving  the  linear  problem  Equa¬ 
tion  (D.3)  to  obtain  the  first  iterates  ai(s)  and  i>i(s).  The  new  Oi(s)  is  then  used  to  solve 
for  02 (s)  and  62(5)  in. 

(D.4) 

etc.  This  method,  due  to  Sanathanan  and  Koerner,  is  often  called  the  SK-iteration.  The  SK- 
iteration  is  not  guaranteed  to  converge.  Whitfield  [48]  has  shown  that  even  if  it  converges, 
it  converges  to  values  that  are  not  the  true  minimum  of  the  NLS  problem.  However,  the 
SK-iteration  is  useful  to  obtain  an  initial  estimate  for  the  minimization  of  Equation  (D.l). 

A  positive  frequency  weighting  functions  w(s)  can  be  included  in  the  cost  function  to 
emphasize  certain  frequency  ranges.  In  addition,  it  is  useful  to  allow  for  known  dynamics 
gk (s)  in  the  measured  transfer  function.  Known  dynamics  include  sensor  and  actuator 
dynamics,  computational  delays,  and  estimates  of  poles  and  zeros  determined  from  Bode 
plots.  The  final  cost  function  thus  has  the  form 

c  =  ~9(s)  •  (D-5) 

a(s ) 

LJ  v  1 

Instead  of  using  this  cost  function,  a  criterion  due  to  Sidman  [44]  can  be  used.  Sidman 
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suggests  fitting  the  log  of  the  transfer  function 


cig  =  E 


log  7u\  ~  log^(s) 

a(s) 


=£ 


m 

_ 

log  j 

a(s) 

\9(s ) 

■  +  * 

arg(^y)  -  arg (g{s)) 


(D-6) 


that  is,  we  axe  minimizing  a  cost  function  consisting  of  the  sum  of  the  ratio  of  the  magni¬ 
tudes  and  phase  difference  of  the  transfer  functions.  Pintelon  [39]  has  shown  that  this  is  a 
consistent  estimator,  that  it  is  not  sensitive  to  plant  model  errors,  (that  is,  the  plant  does 
not  belong  to  the  model  set),  and  that  it  is  robust  to  outliers  in  the  measurements. 


The  log-criterion  cig  is  not  without  problems.  The  phase  is  in  the  range  [—  it,  7r]  and 
thus  “wraps  around,”  resulting  in  discontinuities  that  complicates  the  optimization  step.  A 
way  around  this  is  to  unwrap  the  phase,  but  this,  too,  has  its  own  complications.  This  cost 
function  can  be  extended  to  include  frequency  weighting  and  known  dynamics  as  before.  In 
addition,  one  can  also  include  relative  weighting  between  the  magnitude  and  phase. 


D.2  Parameterization 

If  the  structure  of  the  system  is  known,  e.g.  simple  real  poles/zeros,  complex  poles/zeros, 
delay  elements,  etc.,  the  polynomials  can  be  written  in  factored  form.  For  example,  we  can 
parameterize  the  numerator  polynomial  as 

b(s)  =  '£bjsj  (D.7) 

j=0 

=  fJ(s  +  2j)II(s2 +  2CjfeWfcS  +  u;*).  (D.8) 

i  k 

The  factored  form  has  better  numerical  properties  than  the  polynomial  form  [39].  However, 
if  the  structure  is  not  know  it  is  better  to  use  the  polynomial  representation  because  it  is 
general. 
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Figure  D.l:  Self-similarity  of  standard  polynomials. 

Numerical  problems  axise  because  the  standard  polynomials  s1,  s2,s3, .  •  •  ,sn  become 
self-similar,  see  Adcock  [1].  An  example  of  the  self-similarity  of  the  standard  polynomials 
is  shown  in  Figure  D.l  for  orders  of  practical  interest.  The  higher  degree  polynomials  have 
very  similar  shapes.  To  improve  the  numerical  properties  of  the  standard  polynomials  the 
frequency  range  is  usually  scaled  to  the  range  [0,1],  as  was  done  in  Figure  D.l.  The  self- 
similarity  problem  can  be  solved  by  using  orthogonal  Chebychev  polynomials  instead  of  the 
standard  polynomials.  In  Figure  D.2  we  show  the  corresponding  Chebychev  polynomials. 
The  SK-iteration  typically  converged  2-6  times  faster  for  the  same  problems  when  orthogonal 
polynomials  were  used  instead  of  the  standard  polynomials.  For  discrete  time  systems  it 
is  not  necessary  to  use  Chebychev  polynomials  as  the  complex  exponentials  axe  orthogonal 
on  the  unit  circle. 

Once  we  have  the  numerator  and  denominator  polynomials  we  must  find  a  state  space 
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Figure  D.2:  Chebychev  polynomials. 

representation.  This  is  done  by  computing  the  Markov  parameters  (or  impulse  response) 
from  the  identified  numerator  and  denominator  coefficients,  constructing  Hankel  matrices, 
and  using  the  singular  value  decomposition  to  find  a  minimal  balanced  realization.  This  is 
a  well  known  procedure  and  is  discussed  in  detail  by  Chen  [4].  This  step  only  works  if  the 
system  is  stable.  For  unstable  systems  the  impulse  response  diverges  and  causes  numerical 
problems.  The  procedure  developed  by  Jacques  [23]  was  found  to  work  well  for  unstable 
systems,  even  though  it  is  only  guaranteed  to  work  for  stable  systems. 


D.3  MIMO  System  Identification 

The  MIMO  case  follows  the  same  procedure  as  outlined  above.  All  we  need  is  find  a 
parameterization  for  the  MIMO  transfer  function.  There  are  many  ways  to  do  this,  see 
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Chen  [4].  The  approach  followed  here  is  simply  a  direct  extension  of  the  SISO  description, 
that  is 


G(s) 


B{s) 

a(s) 


(D.9) 


where  B(s)  is  a  matrix  of  polynomials.  Like  the  SISO  case,  the  polynomials  can  also  be 
written  in  factored  form  if  desired.  The  Markov  parameters  and  state  space  matrices  are 
computed  as  before. 


System  identification  is  typically  done  in  several  iterations, 
follows. 


A  typical  iteration  is  as 


1.  Select  a  parameterization  and  solve  the  nonlinear  minimization  Equation  (D.5)  or 
D.6.  This  step  can  be  repeated  with  different  orders  for  the  numerator  and  denomi¬ 
nator  polynomials.  Poles  and  zeros  obtained  from  earlier  iterations  may  be  fixed  (by 
including  them  in  thereby  reducing  the  dimension  of  the  minimization  problem. 

2.  Compute  the  state  space  description  as  discussed  at  the  end  of  Appendix  D.2. 

3.  Compare  transfer  functions  computed  from  the  state  space  description  with  the  mea¬ 
sured  transfer  functions.  Model  reductions  can  be  applied  in  this  step  if  the  dimension 
of  the  state  vector  is  too  high,  else,  repeat  step  1  with  lower  orders  for  the  polynomials. 
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Appendix  E  Parameter  Identification 

The  parameters  in  the  extended  model  can  be  grouped  into  four  sets 


Geometric  parameters 
Pressure  rise  characteristic 
Total  pressure  loss  parameters 
Steady  state  operating  point 


Mr>  Ms>  Ma>  Mi)  ^1>  ^4>  ^G>  ®m) 
dip  dip 

*•  *■  w  s? 

dir  dls 

r,Tf’  d<p'  d(P’ 

<ps,  Its,  hs,  h,  distortion. 


The  axial  measurement  location  xm  has  been  included  in  the  list  of  geometric  param¬ 
eters  although  it  is,  strictly  speaking,  not  part  of  the  model  but  part  of  the  measurement 
relation.  However,  the  sensor  location  has  a  strong  influence  on  the  magnitudes  of  the 
harmonics,  and  therefore  on  the  magnitudes  of  the  transfer  functions,  so  that  it  has  been 
included  here. 

In  Chapter  4  we  have  seen  that  the  dynamics  are  captured  well  by  the  extended  model 
but  there  were  small  differences  in  transfer  function  magnitudes,  and  pole  and  zero  fre¬ 
quencies.  Some  parameters  used  in  the  model  are  first  order  approximations  or  based  on 
prior  experimental  results.  For  example,  the  rotor  and  stator  fluid  inertias  /ir  and  fis  take 
into  account  only  fluid  that  are  in  the  rotors  and  stators  respectively,  while  it  is  reasonable 
to  assume  that  at  least  some  fluid  in  the  gaps  between  the  rotors  and  stators  will  also  be 
accelerated  and  needs  to  be  included  in  these  inertias.  The  proportionality  constant  Tf  used 
in  the  total  pressure  loss  model  was  determined  by  Haynes  [20]  by  fitting  the  individual 
harmonic  transfer  functions  to  experimental  data. 


173 


The  parameters  that  we  would  like  to  determine  are  //r,  /is,  /ia,  xm,  r,  Tf,  || , 


3-0 


The  actuator  fluid  inertia  /za  and  slope  —  have  a  strong  influence  on  the  dominant  zero. 

d'y 

In  modelling  the  unsteady  viscous  effects  it  was  assumed  that  the  losses  divide  between 
the  rotors  and  stators  in  proportion  to  the  reaction  r.  The  time  constants  associated  with 
the  loss  dynamics  were  assumed  to  equal  the  corresponding  flow-through  times  multiplied 
by  the  proportionality  constant  Tf,  see  Equations  (3.61)  and  (3.62).  By  estimating  these 
two  parameters  we  can  verify  our  assumptions.  Even  though  the  slope  —  is  available 
from  the  characteristic,  it  was  found  that  it  is  better  to  estimate  it  as  well.  The  poles 
are  very  sensitive  to  the  slope  and  small  errors  in  it  leads  to  biased  results  for  the  other 
parameters.  The  throttle  constant  is  not  known  and  must  also  be  estimated.  Initially, 
the  parameter  bc  was  also  estimated  and  was  found  to  agree  with  the  geometrical  value  so 
it  is  not  estimated.  We  will  denote  the  list  of  parameters  to  be  estimated  by  the  vector 


c  _  r  dip  dip  ]T 

S  —  LM r?  A^s >  M a?  ^m>  7*,  Tf,  ^ &tj  . 


(E.l) 


Instead  of  using  steady  state  measurements  to  identify  the  required  parameters,  we 
will  determine  them  directly  from  the  transfer  functions.  In  Appendix  C.l  we  saw  that 
the  transfer  functions  can  be  measured  with  high  accuracy  and  thus  we  believe  that  the 
approach  proposed  here  will  result  in  better  estimates  of  the  parameters. 


The  linearized  model  depends  on  the  steady  state  flow  (f>i3)(9 )  which  we  do  not  have,  so 

d'lb 

all  the  quantities,  e.g.,  the  slope  and  steady  state  losses  and  /ss,  are  unknown.  If  we 

d(p 

restrict  ourselves  to  the  uniform  flow  case  these  quantities  are  known  or  can  be  computed 
for  the  specific  operating  point,  so  we  will  use  the  uniform  flow  transfer  functions  to  do  the 
parameter  identification. 


The  same  set  of  parameters  must  fit  all  the  individual  transfer  functions,  not  just  one,  so 
instead  of  using  the  transfer  function  of  an  individual  harmonic  to  identify  the  parameters, 
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we  will  use  all  the  transfer  functions  simultaneously.  Using  all  the  transfer  functions  provide 
us  with  more  independent  data  points  and  thus  the  standard  deviations  of  the  estimates 
are  reduced. 

A  standard  least  squaxes  procedure  is  used.  The  cost  function  c  to  be  minimized  is 

3 

C=  (5nn(*w;£) -p„n(iu>))2  (E.2) 

n=0  i 

where  gnn(iu)] £)  is  the  transfer  function  form  S'jn  to  8(f>n  and  depends  on  the  parameter 
vector  £,  and  gnn{iu)  is  the  measured  transfer  function.  Because  we  are  using  the  uniform 
flow  transfer  functions  only  the  diagonal  elements  of  the  transfer  function  matrix  are  used, 
the  others  are  zero.  The  transfer  functions  gnn(i<jj\  £)  is  highly  nonlinear  in  the  argument 
£  and  a  good  initial  estimate  is  necessary  to  ensure  successful  minimization  of  the  cost 
function.  The  nominal  geometric  and  steady  state  values  were  found  to  be  suitable  for  this 
purpose. 

The  nominal  and  identified  parameters  are  listed  in  the  table  below.  In  this  table  the 
first  row  shows  the  elements  of  the  initial  estimate  f0  and  the  second  row  the  identified 
elements  £. 
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0.68 

0.33 

0.29 

-0.60 

0.75 

1.50 

-0.04 

-0.26 

8.00 

c 

0.58 

0.68 

0.15 

-0.47 

0.76 

1.53 

-0.12 

-0.16 

6.25 

Figure  E.l  shows  the  experimentally  measured  gnn{iuj),  nominal  gnn(iu;  £0),  and  iden¬ 
tified  gnn(iu\ £)  transfer  functions.  All  the  magnitudes  of  the  nominal  transfer  functions 
(dashed  lines)  are  larger  than  the  measured  magnitudes  (dots),  and  the  peaks  (poles)  of  the 
nominal  transfer  function  magnitudes  are  at  slightly  higher  frequencies.  Both  these  differ- 
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ences  are  absent  in  the  identified  transfer  functions  (solid  lines)  and  the  identified  transfer 
functions  clearly  have  a  better  fit. 

The  estimated  rotor  fluid  inertia  /ir  is  15%  smaller  than  the  geometrical  value  while 
the  estimated  stator  fluid  inertia  fis  is  almost  twice  as  large  as  the  geometrical  value.  This 
suggests  that  the  fluid  in  the  gaps  should  be  lumped  with  those  in  the  stators.  If  the  lengths 
of  the  gaps  between  the  rotors  and  stators  and  between  consecutive  stages  are  added  to  that 
of  the  stators,  we  get  an  effective  stator  fluid  inertia  equal  to  0.676.  This  agrees  very  well 
with  the  identified  value  0.68. 

Longley  [25]  used  the  basic  steady  state  model  (no  losses)  to  compute  /zr  and  find  it 
should  be  approximately  twice  the  geometrical  value.  Hynes  modelled  losses  by  modifying 
the  effective  rotor  fluid  inertia.  All  these  are  based  on  steady  state  measurements.  The 
transfer  functions  are  a  direct  measurement  of  the  dynamic  behavior,  the  parameters  of 
which  we  are  trying  to  identify,  so  the  estimates  obtained  from  the  transfer  functions  are 
superior. 

The  actuator  fluid  inertia  /ia  is  smaller  than  the  geometrical  value.  A  possible  expla¬ 
nation  may  be  the  low  solidity  of  the  AG  Vs. 

The  estimated  measurement  location  xm  is  smaller  than  the  geometrical  value  and 
puts  the  source  of  the  perturbations  towards  the  end  of  the  first  stage.  This  agrees  with 
measurements  by  Haynes  [20]  who  found  that  the  different  harmonics  were  strongest  in  the 
first  stage.  The  geometrical  value  was  obtained  by  measuring  the  distance  between  the  inlet 
of  the  first  rotor  and  the  sensors. 

The  reaction  r  and  constant  of  proportionality  Tf  used  in  the  model  of  the  unsteady 
viscous  effects  agree  well  with  the  nominal  values.  Haynes  [20]  estimated  Tf  on  the  same 
compressor  and  his  estimate  was  used  as  the  nominal  value. 
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The  estimated  slope  =  —0.12  corresponds  to  <f>  =  0.473  on  the  constant  speed  line 
o<p 

e  dip  . 

(nominal  <f>  =  0.470)  and  is  within  measurement  error.  The  smaller  value  of  —  indicates 
that  the  actuators  are  less  effective  than  the  model  predicts  and  may  be  a  result  of  the 
simplified  modelling  of  the  AGVs.  In  the  derivation  of  AGV  relations  we  have  assumed  a 
continuum  of  blades  while  there  are  only  12  around  the  annulus. 


We  could  follow  this  parameter  identification  step  with  a  second  step  during  which 
we  use  the  geometric  parameters  that  we  have  just  identified  and  estimate  the  pressure 
rise  characteristic  tp(<p)  over  the  full  range  of  the  distorted  flow  (f)(0)-  However,  to  do  so 
we  need  <p (s3)(0),  the  steady  flow  through  the  compressor,  but  we  only  have  the  upstream 
measurement  ^(O )•  One  solution  is  to  identify  <f>s\0)  as  well.  Estimating  both  ip(4>)  and 
<f>s3)(6)  led  to  divergence  of  the  least  squares  minimization  step.  This  is  probably  because  all 
the  steady  state  and  pressure  rise  parameters  are  now  free  variables  and  the  optimization 
problem  may  not  be  well  posed  any  more.  Therefore,  determination  of  the  characteristic 
over  the  full  range  of  the  nonuniform  flow  is  still  an  open  problem. 
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